mintwo.cpp

最小二乘法

#include "stdafx.h"
#include "math.h"
#include "stdio.h"
#define ABS(x) (x)>0?(x):-(x)
#define SWAP(a,b) {temp=(a);(a)=(b);(b)=temp;}

void solve(double **a,double *b,double *x,int n)
{
	int i,j,k,ik;
	double mik,temp;
	for(k=0;k<n;k++)
	{
		mik=-1;
		for(i=k;i<n;i++)
			if(ABS(a[i][k])>mik)
			{
				mik=ABS(a[i][k]);
				ik=i;
			}
		for(j=k;j<n;j++)
			SWAP(a[ik][j],a[k][j]);
		SWAP(b[k],b[ik]);
		b[k]/=a[k][k];
		for(i=n-1;i>=k;i--)
			a[k][i]/=a[k][k];
		for(i=k+1;i<n;i++)
		{
			b[i]-=a[i][k]*b[k];
			for(j=n-1;j>=k;j--)
				a[i][j]-=a[i][k]*a[k][j];
		}
	}
	for(i=n-1;i>=0;i--)
	{
		x[i]=b[i];
		for(j=i+1;j<n;j++)
			x[i]-=a[i][j]*x[j];
	}
}

void linear(double **x,double *y,double *beta,int n,int p)
{
	double **a,*b;
	int i,j,k;
	a=new double*[p];
	for(i=0;i<p;i++)
		a[i]=new double[p];
	for(i=0;i<p;i++)
		for(j=0;j<p;j++)
		{
			a[i][j]=0;
			for(k=0;k<n;k++)
				a[i][j]+=x[k][i]*x[k][j];
		}
	b=new double[p];
	for(i=0;i<p;i++)
	{
		b[i]=0;
		for(j=0;j<n;j++)
			b[i]+=x[j][i]*y[j];
	}
	solve(a,b,beta,p);
	for(i=0;i<p;i++)
		delete a[i];
	delete a,b;
}

void polyfit(double *x,double *y,double *beta,int n,int p)
{
	int i,j;
	double **xx;
	xx=new double*[n];
	for(i=0;i<n;i++)
	{
		xx[i]=new double[p];
		for(j=0;j<p;j++)
			xx[i][j]=pow(x[i],j);
	}
	linear(xx,y,beta,n,p);
	for(i=0;i<n;i++)
		delete xx[i];
	delete xx;
}

void polyfitw(double *x,double *y,double *w,double *beta,int n,int p)
{
	int i,j;
	double **xx;
	xx=new double*[n];
	for(i=0;i<n;i++)
	{
		w[i]=sqrt(w[i]);
		xx[i]=new double[p];
		for(j=0;j<p;j++)
			xx[i][j]=pow(x[i],j)*w[i];
		y[i]*=w[i];
	}
	linear(xx,y,beta,n,p);
	for(i=0;i<n;i++)
		delete xx[i];
	delete xx;
}

void call1()
{
	FILE *fp;
	int i,n;
	double *x,*y;
	fp=fopen("d:\\data.txt","r");
	fscanf(fp,"%d",&n);
	x=new double[n];
	y=new double[n];
	for(i=0;i<n;i++)
	{
		fscanf(fp,"%lf",x+i);
		fscanf(fp,"%lf",y+i);
	}
	fclose(fp);

	int p=4;
	double *beta;
	beta=new double[p];
	polyfit(x,y,beta,n,p);
	for(i=0;i<p;i++)
		printf("%6.4f  ",beta[i]);
	delete x,y,beta;
}

void call2()
{
	FILE *fp;
	int i,n;
	double *x,*y,*w;
	fp=fopen("d:\\data369.txt","r");
	fscanf(fp,"%d",&n);
	x=new double[n];
	y=new double[n];
	w=new double[n];
	for(i=0;i<n;i++)
	{
		fscanf(fp,"%lf",x+i);
		fscanf(fp,"%lf",y+i);
		fscanf(fp,"%lf",w+i);
	}
	fclose(fp);

	int p=4;
	double *beta;
	beta=new double[p];
	polyfitw(x,y,w,beta,n,p);
	for(i=0;i<p;i++)
		printf("%15.13f  ",beta[i]);
	delete x,y,w,beta;
}

#include "stdlib.h"
void call()
{
	int i,n=1000,p=7;
	double *x,*y,*w,*beta;
	x=new double[n];
	y=new double[n];
	w=new double[n];
	beta=new double[p];
	for(i=0;i<n;i++)
	{
		x[i]=i*1.0/n;
		y[i]=sin(x[i]);//1.0*rand()/RAND_MAX;
		w[i]=fabs(exp(x[i]));
	}
	polyfitw(x,y,w,beta,n,p);
	for(i=0;i<p;i++)
		printf("%10.8f\n",beta[i]);
	delete x,y,w,beta;
}

double fmin(double **a,double *b,double c,double *x,int n)
{
	solve(a,b,x,n);
	double rt=c;
	for(int i=0;i<n;i++)
		rt-=b[i]*x[i];
	return rt;
}

void main()
{
	int i,j,n;
	double **a,*b,c,*x,f;
	FILE *fp;
	fp=fopen("d:\\fmin.txt","r");
	fscanf(fp,"%d",&n);
	a=new double*[n];
	b=new double[n];
	x=new double[n];
	for(i=0;i<n;i++)
		a[i]=new double[n];
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
			fscanf(fp,"%lf",a[i]+j);
		fscanf(fp,"%lf",b+i);
	}
	fscanf(fp,"%lf",&c);
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)
			printf("%6.4f  ",a[i][j]);
		printf("%6.4f\n",b[i]);
	}
	printf("%6.4f\n",c);
	f=fmin(a,b,c,x,n);
	for(i=0;i<n;i++)
		printf("%6.4f  ",x[i]);
	printf("\n%6.4f\n",f);
	fclose(fp);
}