Метод исключения, страница 3

double *Y, *Y_new, *F, **Teta;

double RSS, RSS_new, Sigma2, y_sr;

double *Cm, *F3;

short Iter=0, j_is, **Model, *Preo;

void main () {

short i, j, k, jmax;

double max;

FILE *fx;

fx = fopen( "x.txt", "r" );

fscanf( fx, "%d", &M_Glob );

fscanf( fx, "%d", &n );

M_Glob++;

m = M_Glob;

X = new double* [n];

for ( i=0; i<n; i++ )

X[i] = new double [m];

X_new = new double* [n];

for ( i=0; i<n; i++ )

X_new[i] = new double [m-1];

A = new double* [m];

for ( i=0; i<m; i++ )

A[i] = new double [m];

Teta = new double* [m];

for ( i=0; i<m; i++ )

Teta[i] = new double [m];

Y = new double [n];

Y_new = new double [n];

F = new double [m];

Cm = new double [m];

F3 = new double [m];

Model = new short* [m];

for ( i=0; i<m; i++ )

Model[i] = new short [m];

Preo = new short [m];

for ( i=0; i<n; i++ ){

for ( j=0; j<m-1; j++ )

fscanf( fx, "%lf", &(X[i][j]) );

X[i][m-1]=1;

fscanf( fx, "%lf", &(Y[i]) );

}

fclose( fx );

for ( i=0; i<M_Glob; i++ ) {

Preo[i] = i;

for ( j=0; j<M_Glob; j++ )

Model[i][j] = 1;

}

//////////////////////////////////////////////////////////////////////////

for ( Iter=0; Iter<M_Glob; Iter++ ) {

for ( i=0; i<m; i++ )

for ( j=0; j<m; j++ ) {

A[i][j] = 0;

for ( k=0; k<n; k++ )

A[i][j] += X[k][i]*X[k][j];

}

for ( i=0; i<m; i++ ) {

F[i]=0;

for ( j=0; j<n; j++ )

F[i] += X[j][i]*Y[j];

}

for(i=0;i<m;i++) {

for( jmax=i,max=fabs(A[i][i]),j=i+1;j<m;j++)

if(fabs(A[j][i])>max) {

max=fabs(A[j][i]);

jmax=j;

}

if(jmax>i) {

for(j=i;j<m;j++) {

max=A[i][j];

A[i][j]=A[jmax][j];

A[jmax][j]=max;

}

max=F[i];

F[i]=F[jmax];

F[jmax]=max;

}

F[i]/=A[i][i];

for(j=m-1;j>=i;j--)  A[i][j]/=A[i][i];

for(j=i+1;j<m;j++) {

F[j]-=F[i]*A[j][i];

for(int k=m-1;k>=i;k--) A[j][k]-=A[i][k]*A[j][i];

}

}

for(i=m-1;i>=0;i--) {

for(j=0;j<i;j++) {

F[j]-=A[j][i]*F[i];                                         

}            

}

/////////////////////////////////////////////

for ( j=0; j<m; j++ )

Teta[Iter][j] = F[j];

for ( i=0; i<n; i++ ) {

Y_new[i]=0;

for ( j=0; j<m; j++ )

Y_new[i] += X[i][j]*F[j];

}

RSS = 0;

for ( i=0; i<n; i++ )

RSS += (Y_new[i]-Y[i])*(Y_new[i]-Y[i]);

if ( Iter==0 ) {

Sigma2=RSS/(n-m);

}

Cm[Iter]=RSS/Sigma2+2*m-n;

m--;

for (       j_is=0; j_is<m+1; j_is++ ) {

for ( i=0; i<n; i++ )

for ( j=0; j<m; j++ )

if ( j<j_is )           X_new[i][j]=X[i][j];

else                                     X_new[i][j]=X[i][j+1];

//////////////////////////////////////////////////////////////////////////

for ( i=0; i<m; i++ )

for ( j=0; j<m; j++ ) {

A[i][j] = 0;

for ( k=0; k<n; k++ )

A[i][j] += X_new[k][i]*X_new[k][j];

}

for ( i=0; i<m; i++ ) {

F[i]=0;

for ( j=0; j<n; j++ )

F[i] += X_new[j][i]*Y[j];

}

for(i=0;i<m;i++) {

for( jmax=i,max=fabs(A[i][i]),j=i+1;j<m;j++)

if(fabs(A[j][i])>max) {

max=fabs(A[j][i]);

jmax=j;

}

if(jmax>i) {

for(j=i;j<m;j++) {

max=A[i][j];

A[i][j]=A[jmax][j];

A[jmax][j]=max;

}

max=F[i];

F[i]=F[jmax];

F[jmax]=max;

}

F[i]/=A[i][i];

for(j=m-1;j>=i;j--)  A[i][j]/=A[i][i];

for(j=i+1;j<m;j++) {

F[j]-=F[i]*A[j][i];

for(int k=m-1;k>=i;k--) A[j][k]-=A[i][k]*A[j][i];

}

}

for(i=m-1;i>=0;i--) {

for(j=0;j<i;j++) {

F[j]-=A[j][i]*F[i];                                         

}            

}

/////////////////////////////////////////////

for ( i=0; i<n; i++ ) {

Y_new[i]=0;

for ( j=0; j<m; j++ )

Y_new[i] += X_new[i][j]*F[j];

}

RSS_new = 0;

for ( i=0; i<n; i++ )

RSS_new += (Y_new[i]-Y[i])*(Y_new[i]-Y[i]);

F3[j_is] = (RSS_new-RSS)/RSS_new/(n-(m+1)+1);

}

max = F3[0];

j_is=0;

for ( i=0; i<m+1; i++ )

if ( max>F3[i] ) {

max = F3[i];

j_is = i;

}

for ( i=0; i<n; i++ )

for ( j=j_is; j<m; j++ )

X[i][j]=X[i][j+1];

for ( j=Iter+1; j<M_Glob; j++ )

Model[j][ Preo[j_is] ]=0;

for ( j=j_is; j<m; j++ )

Preo[j] = Preo[j+1];

}

for ( i=0; i<n; i++ )

delete[] X[i];

delete[] X;

for ( i=0; i<n; i++ )

delete[] X_new[i];

delete[] X_new;

for ( i=0; i<M_Glob; i++ )

delete[] A[i];

delete[] A;

delete[] Y;

delete[] Y_new;

delete[] F;

delete[] F3;

max = Cm[0];

j_is=0;

for( i=1; i<M_Glob; i++ )

if ( max>Cm[i] ) {

max = Cm[i];

j_is = i;

}

delete[] Cm;

for ( j=0, i=0; j<M_Glob; j++ )

if ( Model[j_is][j] ) {

if ( j!=0 )

printf( " + " );

if ( j!=M_Glob-1 ) printf( "%5.3lf*X%d", Teta[j_is][i], j+1 );

else printf( "%5.3lf", Teta[j_is][i] );

i++;

}

printf( "\n\n" );

for ( i=0; i<M_Glob; i++ )

delete[] Teta[i];

delete[] Teta;

for ( i=0; i<M_Glob; i++ )

delete[] Model[i];

delete[] Model;

delete[] Preo;

}