Changeset 39246 for trunk/Ohana/src/relastro/src/FitPM_IRLS.c

Timestamp:

Dec 11, 2015, 9:12:57 AM (11 years ago)

Author:

eugene

Message:

adding IRLS versions of fitpm and fitplx

File:

• trunk/Ohana/src/relastro/src/FitPM_IRLS.c (modified) (6 diffs)

trunk/Ohana/src/relastro/src/FitPM_IRLS.c

@@ -11 +11 @@
   int i,j;
 
-  int p   = 4;
-  int n   = 2 * Npoints;
-  int dof = n - p;
-  
-  // data->A,B,Cov,Beta,Beta_prev are allocated outside by FitAstromDataInit()
-  // points->Wx,Wy,rx,ry,u are elements of FitAstromPoint
+  int Ndof = 2 * Npoints - data->Nterms;
   
   // Convert the measurement errors into initial weights.
   for (i = 0; i < Npoints; i++) {
-    // points[i].Wx = 1 / points[i].dX;
-    // points[i].Wy = 1 / points[i].dY;
     points[i].Wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1 / SQ(points[i].dX);
     points[i].Wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1 / SQ(points[i].dY);
@@ -39 +32 @@
     sigma_ols += SQ(points[i].rx) + SQ(points[i].ry);
   }
-  sigma_ols = sqrt(sigma_ols / dof);
+  sigma_ols = sqrt(sigma_ols / (float)Ndof);
 
   // Save OLS covariance and Beta (solution vector, which is actually also saved in fit)
@@ -101 +94 @@
   }
 
-  // this section calculates the formal error on the weighted fit using the covariance values
+  // calculate the weight thresholds to mask the bad points:
   double Sum_Wx = 0.0;
   double Sum_Wy = 0.0;
-  if (1) {
+  for (i = 0; i < Npoints; i++) {
+    points[i].Wx = weight_cauchy(points[i].rx / points[i].dX);
+    points[i].Wy = weight_cauchy(points[i].ry / points[i].dY);
+    
+    Sum_Wx += points[i].Wx;
+    Sum_Wy += points[i].Wy;
+  }
+  double WxThreshold = WEIGHT_THRESHOLD * Sum_Wx / (1.0 * Npoints);
+  double WyThreshold = WEIGHT_THRESHOLD * Sum_Wy / (1.0 * Npoints);
+
+  // set a mask (which can be used by the bootstrap resampling analysis)
+  for (i = 0; i < Npoints; i++) {
+    // keep if either is above threshold?
+    // drop if either is below threshold?
+    // points are marked as keep by default
+    if ((points[i].Wx < WxThreshold) || (points[i].Wy < WyThreshold)) {
+      points[i].mask = 1; // keep point if mask == 0
+    }
+  }
+
+  // this section calculates the formal error on the weighted fit using the covariance values
+  // NOTE EAM: in tests (fitpm.c), they seem to be too large by a factor of ~5.37
+  if (data->getError) {
     double ax = 0.0, ay = 0.0;
     double bx = 0.0, by = 0.0;
-    double lambda = 0.0;
+
     for (i = 0; i < Npoints; i++) {
-      points[i].Wx = weight_cauchy(points[i].rx / points[i].dX);
-      points[i].Wy = weight_cauchy(points[i].ry / points[i].dY);
-    
       ax += dpsi_cauchy(points[i].rx / points[i].dX);
       ay += dpsi_cauchy(points[i].ry / points[i].dY);
@@ -117 +129 @@
       bx += SQ(points[i].Wx);
       by += SQ(points[i].Wy);
-
-      Sum_Wx += points[i].Wx;
-      Sum_Wy += points[i].Wy;
     }
     ax /= 1.0 * Npoints;  // mean(psi_dot(r))
     ay /= 1.0 * Npoints; 
-    bx /= 1.0 * (Npoints - p); // mean(psi^2(r)) * (N / (N-p))
-    by /= 1.0 * (Npoints - p);
-  
-    double sigma_robust_x = lambda * sqrt(bx) * sigma_hat * 2.385 / ax;
-    double sigma_robust_y = lambda * sqrt(by) * sigma_hat * 2.385 / ay;
+    bx /= 1.0 * (Npoints - data->Nterms); // mean(psi^2(r)) * (N / (N-p))
+    by /= 1.0 * (Npoints - data->Nterms);
+  
+    double lambda_x = 1.0 + (data->Nterms / Npoints) * (1 - ax) / ax;
+    double lambda_y = 1.0 + (data->Nterms / Npoints) * (1 - ay) / ay;
+  
+    double sigma_robust_x = lambda_x * sqrt(bx) * sigma_hat * 2.385 / ax;
+    double sigma_robust_y = lambda_y * sqrt(by) * sigma_hat * 2.385 / ay;
 
     // This is actually sigma^2, as that's the factor in the covariance (dumouchel 4.1)
-    double sigma_final_x  = MAX(SQ(sigma_robust_x), (n * SQ(sigma_robust_x) + SQ(p * sigma_ols)) / (n + SQ(p)));
-    double sigma_final_y  = MAX(SQ(sigma_robust_y), (n * SQ(sigma_robust_y) + SQ(p * sigma_ols)) / (n + SQ(p)));
+    double sigma_final_x  = MAX(SQ(sigma_robust_x), (2 * Npoints * SQ(sigma_robust_x) + SQ(data->Nterms * sigma_ols)) / (2 * Npoints + SQ(data->Nterms)));
+    double sigma_final_y  = MAX(SQ(sigma_robust_y), (2 * Npoints * SQ(sigma_robust_y) + SQ(data->Nterms * sigma_ols)) / (2 * Npoints + SQ(data->Nterms)));
 
     fit[0].dRo = sqrt(data->Cov[0][0]);
@@ -144 +156 @@
   }
 
-  // set a mask (which can be used by the bootstrap resampling analysis)
-  double WxThreshold = WEIGHT_THRESHOLD * Sum_Wx / (1.0 * Npoints);
-  double WyThreshold = WEIGHT_THRESHOLD * Sum_Wy / (1.0 * Npoints);
-
-  for (i = 0; i < Npoints; i++) {
-    // keep if either is above threshold?
-    // drop if either is below threshold?
-    // points are marked as keep by default
-    if ((points[i].Wx < WxThreshold) || (points[i].Wy < WyThreshold)) {
-      points[i].mask = 1;
-    }
-  }
-
   // (optionally) add up the chi square for the fit, only counting the unmasked points
   double chisq = 0.0;
@@ -162 +161 @@
   for (i = 0; i < Npoints; i++) {
     if (points[i].mask) continue;
+    fit[0].Nfit ++;
       
-    double Xf = fit[0].Ro + fit[0].uR*points[i].T;
-    double Yf = fit[0].Do + fit[0].uD*points[i].T;
-    double wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dX);
-    double wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dY);
-    chisq += SQ(points[i].X - Xf) * wx;
-    chisq += SQ(points[i].Y - Yf) * wy;
-    fit[0].Nfit ++;
-  }
-    
-  // the reduced chisq is divided by (Ndof = 2*Nfit - 4)
-  fit[0].chisq = chisq / (2.0*fit[0].Nfit - 4.0);
+    if (data->getChisq) {
+      double Xf = fit[0].Ro + fit[0].uR*points[i].T;
+      double Yf = fit[0].Do + fit[0].uD*points[i].T;
+      double wx = (fabs(points[i].dX) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dX);
+      double wy = (fabs(points[i].dY) < 0.0001) ? 1.0 : 1.0 / SQ(points[i].dY);
+      chisq += SQ(points[i].X - Xf) * wx;
+      chisq += SQ(points[i].Y - Yf) * wy;
+    }
+  }
+    
+  // the reduced chisq is divided by (Ndof = 2*Nfit - Nterms)
+  fit[0].chisq = chisq / (2.0*fit[0].Nfit - data->Nterms);
 
   return (TRUE);