/*************************************************** * Program to demonstrate least squares * * polynomial fitting subroutine * * ------------------------------------------------ * * Reference: BASIC Scientific Subroutines, Vol. II * * By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1].* * * * C++ version by J-P Moreau, Paris * * (www.jpmoreau.fr) * * ------------------------------------------------ * * SAMPLE RUN: * * * * LEAST SQUARES POLYNOMIAL FITTING * * * * What is the order of the fit : 4 * * What is the error reduction factor: 0 * * How many data pooints are there : 11 * * * * Input the data points as prompted: * * * * 1 X, Y = 1,1 * * 2 X, Y = 2,2 * * 3 X, Y = 3,3 * * 4 X, Y = 4,4 * * 5 X, Y = 5,5 * * 6 X, Y = 6,6 * * 7 X, Y = 7,7 * * 8 X, Y = 8,8 * * 9 X, Y = 9,9 * * 10 X, Y = 0,0 * * 11 X, Y = 1,1 * * * * Coefficients are: * * 0 -0.000000 * * 1 1.000000 * * 2 -0.000000 * * 3 0.000000 * * 4 -0.000000 * * * * Standard deviation = 0.000000 * * * ***************************************************/ #include #include #define SIZE 25 typedef double TAB[SIZE+1]; int i,l,m,n; double dd,e1,vv; TAB x,y,v,a,b,c,d,c2,e,f; /**************************************************************** * LEAST SQUARES POLYNOMIAL FITTING SUBROUTINE * * ------------------------------------------------------------- * * This program least squares fits a polynomial to input data. * * forsythe orthogonal polynomials are used in the fitting. * * The number of data points is n. * * The data is input to the subroutine in x[i], y[i] pairs. * * The coefficients are returned in c[i], * * the smoothed data is returned in v[i], * * the order of the fit is specified by m. * * The standard deviation of the fit is returned in d. * * There are two options available by use of the parameter e: * * 1. if e = 0, the fit is to order m, * * 2. if e > 0, the order of fit increases towards m, but will * * stop if the relative standard deviation does not decrease * * by more than e between successive fits. * * The order of the fit then obtained is l. * ****************************************************************/ void LS_POLY() { //Labels: e10,e15,e20,e30,e50,fin; int i,l2,n1; double a1,a2,b1,b2,c1,d1,f1,f2,v1,v2,w; n1 = m + 1; v1 = 1e7; // Initialize the arrays for (i = 1; i < n1+1; i++) { a[i] = 0; b[i] = 0; f[i] = 0; }; for (i = 1; i < n+1; i++) { v[i] = 0; d[i] = 0; } d1 = sqrt(n); w = d1; for (i = 1; i < n+1; i++) { e[i] = 1 / w; } f1 = d1; a1 = 0; for (i = 1; i < n+1; i++) { a1 = a1 + x[i] * e[i] * e[i]; } c1 = 0; for (i = 1; i < n+1; i++) { c1 = c1 + y[i] * e[i]; } b[1] = 1 / f1; f[1] = b[1] * c1; for (i = 1; i < n+1; i++) { v[i] = v[i] + e[i] * c1; } m = 1; e10: // Save latest results for (i = 1; i < l+1; i++) c2[i] = c[i]; l2 = l; v2 = v1; f2 = f1; a2 = a1; f1 = 0; for (i = 1; i < n+1; i++) { b1 = e[i]; e[i] = (x[i] - a2) * e[i] - f2 * d[i]; d[i] = b1; f1 = f1 + e[i] * e[i]; } f1 = sqrt(f1); for (i = 1; i < n+1; i++) e[i] = e[i] / f1; a1 = 0; for (i = 1; i < n+1; i++) a1 = a1 + x[i] * e[i] * e[i]; c1 = 0; for (i = 1; i < n+1; i++) c1 = c1 + e[i] * y[i]; m = m + 1; i = 0; e15: l = m - i; b2 = b[l]; d1 = 0; if (l > 1) d1 = b[l - 1]; d1 = d1 - a2 * b[l] - f2 * a[l]; b[l] = d1 / f1; a[l] = b2; i = i + 1; if (i != m) goto e15; for (i = 1; i < n+1; i++) v[i] = v[i] + e[i] * c1; for (i = 1; i < n1+1; i++) { f[i] = f[i] + b[i] * c1; c[i] = f[i]; } vv = 0; for (i = 1; i < n+1; i++) vv = vv + (v[i] - y[i]) * (v[i] - y[i]); //Note the division is by the number of degrees of freedom vv = sqrt(vv / (n - l - 1)); l = m; if (e1 == 0) goto e20; //Test for minimal improvement if (fabs(v1 - vv) / vv < e1) goto e50; //if error is larger, quit if (e1 * vv > e1 * v1) goto e50; v1 = vv; e20: if (m == n1) goto e30; goto e10; e30: //Shift the c[i] down, so c(0) is the constant term for (i = 1; i < l+1; i++) c[i - 1] = c[i]; c[l] = 0; //l is the order of the polynomial fitted l = l - 1; dd = vv; goto fin; e50: // Aborted sequence, recover last values l = l2; vv = v2; for (i = 1; i < l+1; i++) c[i] = c2[i]; goto e30; fin: ; } void main() { printf(" LEAST SQUARES POLYNOMIAL FITTING\n"); printf("\n What is the order of the fit : "); scanf("%d",&m); printf("\n What is the error reduction factor: "); scanf("%lf",&e1); printf("\n How many data pooints are there : "); scanf("%d",&n); printf("\n\n Input the data points as prompted:\n"); for (i = 1; i < n+1; i++) { printf(" %d X Y = ",i); scanf("%lf %lf",&x[i], &y[i]); } printf("\n"); LS_POLY(); // CALL LEAST SQUARES POLYNOMIAL FITTING SUBROUTINE printf("\n Coefficients are:\n"); for (i = 0; i < l+1; i++) { printf(" %d %f\n",i,c[i]); } printf("\n\n Standard deviation: %f\n\n",dd); } // End of file lsqply.cpp