?? 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 <stdio.h>
#include <math.h>
#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<n+1; i++)
for (j=1; j<m+1; j++)
B[i][j] = A[j][i];
}
// Matrix save (A in C)
void A_IN_C(int n1,int n2) {
int i1,i2;
if (n1 * n2 != 0)
for (i1 = 1; i1<n1+1; i1++)
for (i2 = 1; i2<n2+1; i2++)
C[i1][i2] = A[i1][i2];
}
// Matrix save (B in A)
void B_IN_A(int n1,int n2) {
int i1,i2;
if (n1 * n2 != 0)
for (i1 = 1; i1<n1+1; i1++)
for (i2 = 1; i2<n2+1; i2++)
A[i1][i2] = B[i1][i2];
}
// Matrix save (C in B)
void C_IN_B(int n1,int n2) {
int i1,i2;
if (n1 * n2 != 0)
for (i1 = 1; i1<n1+1; i1++)
for (i2 = 1; i2<n2+1; i2++)
B[i1][i2] = C[i1][i2];
}
// Matrix save (C in A)
void C_IN_A(int n1,int n2) {
int i1,i2;
if (n1 * n2 != 0)
for (i1 = 1; i1<n1+1; i1++)
for (i2 = 1; i2<n2+1; i2++)
A[i1][i2] = C[i1][i2];
}
// Matrix multiplication C = A * B
// A is m1 by n1, B is m2 by n2 ( m2 = n1 )
// Result C is m1 by n2.
void Matmult(int m1,int n1,int n2) {
int i,j,k;
for (i = 1; i < m1+1; i++)
for (j = 1; j < n2+1; j++) {
C[i][j] = 0;
for (k = 1; k < n1+1; k++)
C[i][j] += A[i][k] * B[k][j];
}
}
// Matrix inversion: B = Inv(A) by Gauss-Jordan method
// A and B are n by n matrices
void Matinv() {
// Labels: E10, E20, E30;
int i,j,k; double bb;
for (i = 1; i < n+1; i++) {
for (j = 1; j < n+1; j++) {
B[i][j + n] = 0;
B[i][j] = A[i][j];
}
B[i][i + n] = 1;
}
for (k = 1; k < n+1; k++) {
if (k == n) goto E10;
m = k;
for (i = k + 1; i < n+1; i++)
if (fabs(B[i][k]) > 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<n+2) m=n+2;
printf("\n\n Input the data in (x,y) pairs as prompted\n");
for (i = 1; i < m+1; i++) {
printf(" %d X Y = ",i); scanf("%lf %lf",&X[i],&Y[i]);
}
printf("\n");
// Call coefficients generation subroutine
Gene_Coeffs();
// Call regression subroutine
n++;
LstSqrN();
printf(" The calculated coefficients are:\n");
for (i = 1; i < n+1; i++) {
printf(" %d %f\n",i-1,D[i]);
}
n--;
// Call standard deviation subroutine
Standard_deviation();
printf("\n\n Standard deviation: %f\n\n",dd);
} // end of main()
// End leastsqr.cpp
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -