1 /*
2 * PolynomialFitter.java
3 *
4 * Created on March 20, 2006, 3:40 PM
5 *
6 *
7 */
8
9 package org.lcsim.fit.polynomial;
10
11 import Jama.Matrix;
12 import Jama.QRDecomposition;
13
14 /**
15 * A least-squares fitter for fitting polynomials to 2-D points.
16 * @author Norman Graf
17 */
18 public class PolynomialFitter
19 {
20
21 private PolynomialFit _fit;
22
23 /** Creates a new instance of PolynomialFitter */
24 public PolynomialFitter()
25 {
26 }
27
28 /**
29 * Perform the least-squares fit.
30 * @param x An array of the independent variable, assumed to have no uncertainty.
31 * @param y The array of dependent variables, assumed to depend on x via a polynomial.
32 * @param sigma_y The array of uncertainties in the dependent variable.
33 * @param NP The number of points to fit.
34 * @param order The order of the polynomial to fit to the data.
35 * @return True if the fit is successful.
36 */
37 public boolean fit(double[] x, double[] y, double[] sigma_y, int NP, int order)
38 {
39 double [][] alpha = new double[order][order];
40 double [] beta = new double[order];
41 double term = 0;
42
43 for (int k=0; k < order; k++)
44 {
45 // Only need to calculate the diagonal and upper half
46 // of symmetric matrix.
47 for (int j=k; j < order; j++)
48 {
49 // Calculate the matrix terms over the data points
50 term = 0.0;
51 alpha[k][j] = 0.0;
52 for (int i=0; i < NP; i++)
53 {
54
55 double prod1 = 1.0;
56 // Calculate x^k
57 if ( k > 0) for (int m=0; m < k; m++) prod1 *= x[i];
58
59 double prod2 = 1.0;
60 // Calculate x^j
61 if ( j > 0) for (int m=0; m < j; m++) prod2 *= x[i];
62
63 // Calculate x^k * x^j
64 term = (prod1*prod2);
65
66 if (sigma_y != null && sigma_y[i] != 0.0) term /= (sigma_y[i]*sigma_y[i]);
67
68 alpha[k][j] += term;
69 }
70 alpha[j][k] = alpha[k][j];// C will need to be inverted.
71 }
72
73 for (int i=0; i < NP; i++)
74 {
75 double prod1 = 1.0;
76 if (k > 0) for ( int m=0; m < k; m++) prod1 *= x[i];
77 term = (y[i] * prod1);
78 if (sigma_y != null && sigma_y[i] != 0.0)
79 term /= (sigma_y[i]*sigma_y[i]);
80 beta[k] +=term;
81 }
82 }
83
84 // Use the Jama QR Decomposition classes to solve for
85 // the parameters.
86 Matrix alpha_matrix = new Matrix(alpha);
87 QRDecomposition alpha_QRD = new QRDecomposition(alpha_matrix);
88 Matrix beta_matrix = new Matrix(beta,order);
89 Matrix parameters;
90 try
91 {
92 parameters = alpha_QRD.solve(beta_matrix);
93 }
94 catch (Exception e)
95 {
96 System.out.println("QRD solve failed: "+ e);
97 return false;
98 }
99
100 // The inverse provides the covariance matrix.
101 Matrix covariance = alpha_matrix.inverse();
102
103 _fit = new PolynomialFit(parameters, covariance);
104 return true;
105 }
106
107 /**
108 * Return the result of the fit.
109 * @return A PolynomialFit object representing the result of the fit.
110 */
111 public PolynomialFit getFit()
112 {
113 return _fit;
114 }
115 }