1 /*
2 * SlopeInterceptLineFitter.java
3 *
4 * Created on March 27, 2006, 3:19 PM
5 *
6 * $Id: SlopeInterceptLineFitter.java,v 1.1 2006/03/28 04:51:03 ngraf Exp $
7 */
8
9 package org.lcsim.fit.line;
10 import static java.lang.Math.sqrt;
11 /**
12 * A least-squares fit to a 2-D line in slope-intercept form y=a+bx.
13 * It assumes that there is no error on the x-component and that all
14 * the error can be projected onto the y-component.
15 * @author Norman Graf
16 */
17 public class SlopeInterceptLineFitter
18 {
19 private SlopeInterceptLineFit _fit;
20 /** Creates a new instance of SlopeInterceptLineFitter */
21 public SlopeInterceptLineFitter()
22 {
23
24 }
25
26 /**
27 * Perform the fit.
28 * @param x The array of independent variables.
29 * @param y The array of dependent variables.
30 * @param sigma_y The array of uncertainties on the dependent variables.
31 * @param n The number of points to fit
32 * @return true if the fit is successful.
33 */
34 public boolean fit(double[] x, double[] y, double[] sigma_y, int n)
35 {
36 double sum,sx,sy,sxx,sxy,syy,det;
37 double chisq = 999999.0;
38 int i;
39
40 if (n < 2)
41 {
42 return false; //too few points, abort
43 }
44
45 //initialization
46 sum = sx = sy = sxx = sxy = syy = 0.;
47
48 //find sum , sumx ,sumy, sumxx, sumxy
49 double[] w = new double[n];
50 for (i=0; i<n; ++i)
51 {
52 w[i] = 1/(sigma_y[i]*sigma_y[i]);
53 sum += w[i];
54 sx += w[i]*x[i];
55 sy += w[i]*y[i];
56 sxx += w[i]*x[i]*x[i];
57 sxy += w[i]*x[i]*y[i];
58 syy += w[i]*y[i]*y[i];
59 }
60
61 det = sum*sxx-sx*sx;
62 if (Math.abs(det) < 1.0e-20) return false; //Zero determinant, abort
63
64 //compute the best fitted parameters A,B
65
66 double slope = (sum*sxy-sx*sy)/det;
67 double intercept = (sy*sxx-sxy*sx)/det;
68
69 //calculate chisq-square
70
71 chisq = 0.0;
72 for (i=0; i<n; ++i)
73 {
74 chisq += w[i]*((y[i])-slope*(x[i])-intercept)*((y[i])-slope*(x[i])-intercept);
75 }
76
77 double slopeErr = sqrt(sum/det);
78 double interceptErr = sqrt(sxx/det);
79 double sigab = -sx/det;
80
81 _fit = new SlopeInterceptLineFit(slope, intercept, slopeErr, interceptErr, sigab, chisq, n-2 );
82 return true;
83 }
84
85 /**
86 * Return the fit.
87 * @return The fit results. Returns null if the fit is not successful.
88 */
89 public SlopeInterceptLineFit getFit()
90 {
91 return _fit;
92 }
93 }