/****************************************************** * PROGRAM TO DEMONSTRATE ONE DIMENSIONAL OPERATION * * OF THE MULTI-NONLINEAR REGRESSION SUBROUTINE * * --------------------------------------------------- * * Ref.: BASIC Scientific Subroutines Vol. II, * * By F.R. Ruckdeschel, Byte/McGRAW-HILL,1981 [1]. * * * * C++ version by J-P Moreau, Paris * * (www.jpmoreau.fr) * * --------------------------------------------------- * * SAMPLE RUN: * * * * MULTI-NON LINEAR REGRESSION * * * * Number of data points (mini 3, maxi 25)........: 11 * * Degree of the polynomial to be fitted (maxi 10): 4 * * * * Input the data in (x,y) pairs 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 * * * * The calculated coefficients are: * * 0 -0.000000 * * 1 1.000000 * * 2 -0.000000 * * 3 0.000000 * * 4 -0.000000 * * * * Standard deviation: 0.000000 * * * ******************************************************/ #include #include #define MMAX 25 #define NMAX 10 double D[NMAX+2], X[MMAX+1], Y[MMAX+1]; int i,m,n; double dd,dx,xmn,xmx,xx,ymn,ymx,yy; double Z[MMAX+1][NMAX+2]; // maxi m=25 data points double A[MMAX+1][MMAX+1]; // maxi n=10 order of regression double B[MMAX+1][2*MMAX+1]; double C[MMAX+1][MMAX+1]; // Matrix transpose B = Tr(A) // A and B are n by m matrices void Transpose() { int i,j; for (i=1; i fabs(B[m][k])) m = i; if (m == k) goto E10; for (j = k; j < 2*n+1; j++) { bb = B[k][j]; B[k][j] = B[m][j]; B[m][j] = bb; } E10:for (j = k + 1; j < 2*n+1; j++) B[k][j] /= B[k][k]; if (k == 1) goto E20; for (i = 1; i < k; i++) for (j = k + 1; j < 2*n+1; j++) B[i][j] -= B[i][k] * B[k][j]; if (k == n) goto E30; E20:for (i = k + 1; i < n+1; i++) for (j = k + 1; j < 2*n+1; j++) B[i][j] -= B[i][k] * B[k][j]; } // k loop E30:for (i = 1; i < n+1; i++) for (j = 1; j < n+1; j++) B[i][j] = B[i][j + n]; } // Matinv() // Coefficient matrix generation void Gene_Coeffs() { int i,j; double b; for (i = 1; i < m+1; i++) { b = 1; for (j = 1; j < n + 2; j++) { Z[i][j] = b; b = b * X[i]; } } } // Matinv() // Least squares multidimensional fitting subroutine void LstSqrN() { int i,j,m4,n4; m4 = m; n4 = n; for (i = 1; i < m+1; i++) for (j = 1; j < n+1; j++) A[i][j] = Z[i][j]; Transpose(); // B=Transpose(A) A_IN_C(m,n); // move A to C B_IN_A(n,m); // move B to A C_IN_B(m,n); // move C to B Matmult(n,m,n); // multiply A and B C_IN_A(n,n); // move C to A Matinv(); // B=Inverse(A) m = m4; // restore m B_IN_A(n,n); // move B to A for (i = 1; i < m+1; i++) for (j = 1; j < n+1; j++) B[j][i] = Z[i][j]; Matmult(n,n,m); // multiply A and B C_IN_A(n,m); // move C to A for (i = 1; i < m+1; i++) B[i][1] = Y[i]; Matmult(n,m,1); // multiply A and B // Product C is n by 1 - Regression coefficients are in C(i,1) // and put for convenience into vector D(i). for (i = 1; i < n+1; i++) D[i] = C[i][1]; } // LstSqrN() // Standard deviation calculation for a polynomial fit void Standard_deviation() { int i, j; double b, yy; dd = 0; for (i = 1; i < m+1; i++) { yy = 0; b = 1; for (j = 1; j < n+2; j++) { yy = yy + D[j] * b; b = b * X[i]; } dd += (yy - Y[i]) * (yy - Y[i]); } if (m - n - 1 > 0) dd /= m - n - 1; else dd = 0; dd = sqrt(dd); } void main() { for (i=0; i < NMAX+2; i++) D[i]=0; printf(" MULTI-NONLINEAR REGRESSION\n\n"); printf(" Number of data points (mini 3, maxi 25)........: "); scanf("%d",&m); if (m<3) m=3; if (m>25) m=25; printf(" Degree of the polynomial to be fitted (maxi 10): "); scanf("%d",&n); if (n<1) n=1; if (n>10) n=10; if (m