/****************************************************** * Program to demonstrate the Parafit subroutine * * --------------------------------------------------- * * Reference: BASIC Scientific Subroutines, Vol. II * * By F.R. Ruckdeschel, BYTE/McGRAWW-HILL, 1981 [1]. * * * * C++ Release 1.0 By J-P Moreau, Paris. * * (www.jpmoreau.fr) * * --------------------------------------------------- * * SAMPLE RUN: * * * * Parametric least squares fit * * * * The input data are: * * * * X( 1) = 1.000000 Y( 1) = 0.033702 * * X( 2) = 2.000000 Y( 2) = 0.249029 * * X( 3) = 3.000000 Y( 3) = 0.944733 * * X( 4) = 4.000000 Y( 4) = 1.840089 * * X( 5) = 5.000000 Y( 5) = 1.840089 * * X( 6) = 6.000000 Y( 6) = 0.944733 * * X( 7) = 7.000000 Y( 7) = 0.249029 * * X( 8) = 8.000000 Y( 8) = 0.033702 * * X( 9) = 9.000000 Y( 9) = 0.002342 * * X(10) = 10.000000 Y(10) = 0.000084 * * * * The coefficients are: * * 2.000000 * * 4.500000 * * 3.000000 * * * * The standard deviation of the fit is 0.000000 * * * * The number of iterations was 49 * ******************************************************/ #include #include #define SIZE 25 double A[SIZE],E1[SIZE],X[SIZE],Y[SIZE]; double a0,d,e,ee1,l1,l2,m0,m1,xx,yy; int i,l,m,n; void S500() { //Function subroutine yy = A[1] * exp(-(xx - A[2]) * (xx - A[2]) / A[3]); } // Residual generation subroutine void S200() { int j; l2 = 0; for (j = 1; j < n+1; j++) { xx = X[j]; // Obtain function S500(); l2 = l2 + (Y[j] - yy) * (Y[j] - yy); } d = sqrt(l2 / (n - l)); } /************************************************************** * Parametric least squares curve fit subroutine * * ----------------------------------------------------------- * * This program least squares fits a function to a set of data * * values by successively reducing the variance. Convergence * * depends on the initial values and is not assured. * * n pairs of data values, X(i), Y(i), are given. There are l * * parameters, A(j), to be optimized across. * * Required are initial values for the A(l) and e. Another * * important parameter which affects stability is e1, which is * * initially converted to e1(l) for the first intervals. * * The parameters are multiplied by (1 - e1(i)) on each pass. * **************************************************************/ void Param_LS() { // Labels: e50,e100, fin int i; for (i = 1; i < l+1; i++) E1[i] = ee1; m = 0; //Set up test residual l1 = 1e6; //Make sweep through all parameters e50: for (i = 1; i < l+1; i++) { a0 = A[i]; //Get value of residual A[i] = a0; e100: S200(); //Store result in m0 m0 = l2; //Repeat for m1 A[i] = a0 * (1.0 - E1[i]); S200(); m1 = l2; //Change interval size if called for //If variance was increased, halve E1(i) if (m1 > m0) E1[i] = -E1[i] / 2.0; //If variance was reduced, increase step size by increasing E1(i) if (m1 < m0) E1[i] = 1.2 * E1[i]; //If variance was increased, try to reduce it if (m1 > m0) A[i] = a0; if (m1 > m0) goto e100; } // i loop //End of a complete pass //Test for convergence m++; if (l2==0) goto fin; if (fabs((l1 - l2) / l2) < e) goto fin; //If this point is reached, another pass is called for l1 = l2; goto e50; fin: ; } // Param_LS() void main() { n = 10; l = 3; printf(" Parametric least squares fit\n\n"); printf(" The input data are:\n\n"); for (i = 1; i < n+1; i++) { X[i] = (double) i; Y[i] = 2.0 * exp(-(X[i] - 4.5) * (X[i] - 4.5) / 3.0); printf(" X(%2d) = %9.6f",i,X[i]); printf(" Y(%2d) = %9.6f\n",i,Y[i]); } printf("\n"); e = 0.1; ee1 = 0.5; A[1] = 10.0; A[2] = 10.0; A[3] = 10.0; Param_LS(); printf(" The coefficients are:\n"); printf(" %f\n",A[1]); printf(" %f\n",A[2]); printf(" %f\n",A[3]); printf("\n The standard deviation of the fit is %f\n\n",d); printf(" The number of iterations was %d\n\n",m); } // End parafit.cpp