www.pudn.com > OptLstSqr.rar > lstsqr.cpp
/******************************************************
* 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 *
* *
* C++ version by J-P Moreau, Paris *
* --------------------------------------------------- *
* 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