least squares fitting algorithm.txt

最小二乘法算法源码资料小合集。多种语言编写。

*************************************************************** 
* 本算法用最小二乘法依据指定的M个基函数及N个已知数据进行曲线拟和 
* 输入: m--已知数据点的个数M 
* f--M维基函数向量 
* n--已知数据点的个数N-1 
* x--已知数据点第一坐标的N维列向量 
* y--已知数据点第二坐标的N维列向量 
* a--无用 
* 输出: 函数返回值为曲线拟和的均方误差 
* a为用基函数进行曲线拟和的系数, 
* 即a[0]f[0]+a[1]f[1]+...+a[M]f[M]. 
****************************************************************/ 
double mini_product(int m,double(*f[M])(double),int n,double x[N], 
double y[N],double a[M]) 
{ 
double e,ff,b[M][M],c[M][1]; 
int i,j,k; 

for(j=0;j<m;j++) /*计算最小均方逼近矩阵及常向量*/ 
{ 
for(k=0;k<m;k++) 
{ 
b[j][k]=0.0; 
for(i=0;i<n;i++) 
b[j][k]+=(*f[j])(x)*(*f[k])(x); 
} 
c[j][0]=0.0; 
for(i=0;i<n;i++) 
c[j][0]+=(*f[j])(x)*y; 
} 
gaussian_elimination(m,b,1,c); /*求拟和系数*/ 
for(i=0;i<m;i++) 
a=c[0]; 
e=0.0; 
for(i=0;i<n;i++) /*计算均方误差*/ 
{ 
ff=0.0; 
for(j=0;j<m;j++) 
ff+=a[j]*(*f[j])(x); 
e+=(y-ff)*(y-ff); 
} 
return(e); 
} 




/************************************************************************* 
* 高斯列主元素消去法求解矩阵方程AX=B,其中A是N*N的矩阵,B是N*M矩阵 
* 输入: n----方阵A的行数 
* a----矩阵A 
* m----矩阵B的列数 
* b----矩阵B 
* 输出: det----矩阵A的行列式值 
* a----A消元后的上三角矩阵 
* b----矩阵方程的解X 
**************************************************************************/ 
double gaussian_elimination(int n,double a[M][M],int m,double b[M][1]) 
{ 
int i,j,k,mk; 
double det,mm,f; 
det = 1.0; 
for(k = 0;k<n-1;k++) /*选主元并消元*/ 
{ 
mm=a[k][k]; 
mk = k; 
for(i=k+1;i<n;i++) /*选择第K列主元素*/ 
{ 
if(fabs(mm)<fabs(a[k])) 
{ 
mm = a[k]; 
mk = i; 
} 
} 
if(fabs(mm)<EPS) 
return(0); 
if(mk!=k) /* 将第K列主元素换行到对角线上*/ 
{ 
for(j=k;j<n;j++) 
{ 
f = a[k][j]; 
a[k][j]=a[mk][j]; 
a[mk][j]=f; 
} 
for(j=0;j<m;j++) 
{ 
f = b[k][j]; 
b[k][j]=b[mk][j]; 
b[mk][j]=f; 
} 
det = -det; 
} 
for(i=k+1;i<n;i++) /*将第K列对角线以下消元为零*/ 
{ 
mm = a[k]/a[k][k]; 
a[k]=0.0; 
for(j=k+1;j<n;j++) 
a[j]=a[j]-mm*a[k][j]; 
for(j=0;j<m;j++) 
b[j]=b[j]-mm*b[k][j]; 
} 
det = det*a[k][k]; 
} 
if(fabs(a[k][k])<EPS) 
return 0; 
det=det*a[k][k]; 
for(i=0;i<m;i++) /*回代求解*/ 
{ 
b[n-1]=b[n-1]/a[n-1][n-1]; 
for(j=n-2;j>=0;j--) 
{ 
for(k=j+1;k<n;k++) 
b[j]=b[j]-a[j][k]*b[k]; 
b[j]=b[j]/a[j][j]; 
} 
} 
return(det); 
} 
 
hua1081979 (2006-01-18 20:46:06)  
我用c#实现过，写了很多年了。。。贴上来给你吧。。。
好像是这个样子的。。。
using System;
using System.Drawing;


namespace Eys.CurvePic
{
/// <summary>
/// 该类用于生成图片
/// </summary>
public class MakePic
{
private Bitmap image;
private Graphics g;
private string XValue;
private int x,y;
//private ArrarList XValue=new ArrayList();
private float max;//Y轴刻度的最大值
private int ScaleCount;//X轴刻度的条数

private float ySValue;//Y轴上单个刻度所对应的统计的数值的部份

public MakePic(int x,int y,string xValue)
{
this.x=x;
this.y=y;
this.image=new Bitmap(x+30,y+30);
this.g=Graphics.FromImage(this.image);
this.XValue=xValue;

SetCoordinate();//设置坐标系
}

public void SetImgSave(string filePath)
{
this.image.Save(filePath,System.Drawing.Imaging.ImageFormat.Gif);
}

//传入Y的值，根据Y的值生成曲线
public void SetCurve(float [] YValue,Pen p) //画出曲线
{
int i;
this.max=YValue[0];
for(i=0;i<YValue.Length;i++)
if (YValue[i]>max)
this.max=YValue[i];//求出数组YValue中的最大值

//将实际值与刻值进行换算

ySValue = max/15;//Y轴上单个刻度所对应的统计的数值的部份

int ScaleX=YValue.Length;
this.ScaleCount=ScaleX; 

//处理X轴刻度
int XDate=this.x-20;//实际画出的X轴长
XDate-=20;//留下20的长度不作处理
int SX=XDate/(ScaleX);//求出X轴单个刻度的长
int StartX=30;//开始时的X轴长

//处理Y轴的刻度
int YRen=this.y-20;//实际Y轴长度
YRen-=20;//留下20的长度不作处理
int SY=YRen/15;//Y轴固定画15个刻度，每个刻实际长度
int StartY=this.y-10;//开始时Y轴的开始点

//开始画曲线
//Pen p=new Pen(Color.Red,1.0f);

int EndX = 0;
int EndY = 0;
int CurveX = StartX - SX;
Point [] ppp = new Point[ScaleX];
for(i = 0;i < ScaleX;i++)
{
CurveX += SX;//开点X坐标
//计Y点的开始位置
float ScaleYValue = YValue[i]/ySValue;//求对应Y轴的标筌刻度
int ScaleY = (int)(ScaleYValue*SY);//求对应Y轴的实际长度
int CurveY = StartY-ScaleY;//得出Y轴的实际刻度
if (i == 0)//第一次执行时不画点，只对其进行赋开始点位置
{
EndX = CurveX;
EndY = CurveY;
continue;
}
Point P1 = new Point(EndX,EndY);//定义开始点

Point P2 = new Point(CurveX,CurveY);//定义结束点
//this.g.DrawLine(p,P1,P2);//画直线
ppp[i-1] = P1;
//ppp[i] = P2;

EndX = CurveX;//将开始的点作为下次的起始点
EndY = CurveY;
}
Point P3 = new Point(EndX,EndY);
ppp[i-1] = P3;
this.g.DrawCurve(p,ppp,1);
p.Dispose();
}

/// <summary>
/// 设置坐标系及生成的图像大小
/// </summary>
private void SetCoordinate()
{
//Graphics g=Graphics.FromImage(this.image);
this.g.Clear(Color.White);
Pen p=new Pen(Color.Black,1.0f);

//设置x坐标开始点
Point px1=new Point(0,y-10);
Point px2=new Point(x+10,y-10);
//设置y坐标开始点
Point py1=new Point(30,10);
Point py2=new Point(30,y+10);

this.g.DrawLine(p,px1,px2);//画x坐标
this.g.DrawLine(p,py1,py2);//画y坐标

this.g.DrawString("日期",new Font("Courier New", 10),new SolidBrush(Color.Red),x-50,y+10);//标出X
this.g.DrawString("人",new Font("Courier New", 10),new SolidBrush(Color.Red),0,10);//标出Y
this.g.DrawString("数",new Font("Courier New", 10),new SolidBrush(Color.Red),0,25);
//this.g.DrawString("(0,"+XValue+")",new Font("Courier New", 10),new SolidBrush(Color.Red),30,y);//标出原点

Point XU=new Point(x+5,y-13);//标出X轴的箭头
Point XD=new Point(x+5,y-7);
this.g.DrawLine(p,px2,XU);
this.g.DrawLine(p,px2,XD);

Point YL=new Point(27,15);//标出Y轴的箭头
Point YR=new Point(33,15);
this.g.DrawLine(p,py1,YL);
this.g.DrawLine(p,py1,YR); 
}
/// <summary>
/// 设置坐标系的刻度
/// </summary>
/// <param name="max">Y轴的最大值</param>
/// <param name="ScaleX">X轴刻度的个数</param>
public void SetScale()
{
int ScaleX = this.ScaleCount;
Pen p = new Pen(Color.Blue,1.0f);
int XDate = this.x-20;//实际画出的X轴长
XDate -= 20;//留下20的长度不作处理
int SX = XDate/ScaleX;//求出X轴单个刻度的长
int i;
int YU = this.y-13;
int YD = this.y-7;
int StartX = 30;//开始时的X轴长

Pen pb = new Pen(Color.Gray,1.0f);
//开始标了出X轴的刻度
for(i = 0;i < ScaleX;i++)
{
StartX += SX;
Point Xiu = new Point(StartX,YU);//设X轴的点
Point Xid = new Point(StartX,YD);//设Y轴的点
Point YY = new Point(StartX,30);
this.g.DrawLine(p,Xiu,Xid);
this.g.DrawLine(pb,Xiu,YY);
}

//处理Y轴的刻度
int YRen = this.y-20;//实际Y轴长度
YRen -= 20;//留下20的长度不作处理
int SY = YRen/15;//Y轴固定画15个刻度，每个刻实际长度
int XL = 27;
int XR = 33;
int StartY = y-10;//开始时Y轴的开始点
int ScaleValueY = 0;//Y轴刻度值

//开始画出Y轴刻度
this.g.DrawString("0",new Font("Courier New", 10),new SolidBrush(Color.Blue),0,StartY);//标出0刻度值
for(i = 0;i < 15;i++)
{
StartY -= SY;
Point Yil = new Point(XL,StartY);//设Y轴的点
Point Yir = new Point(XR,StartY);//设Y轴的点
Point XX = new Point(x+10,StartY);
this.g.DrawLine(p,Yil,Yir); 
this.g.DrawLine(pb,Yir,XX);
ScaleValueY = (int)((i+1)*ySValue);
this.g.DrawString(ScaleValueY.ToString(),new Font("Courier New", 10),new SolidBrush(Color.Blue),0,StartY);//标出刻度值
}
}


//据最小二承法生成曲线图
private void LeaseTowRideCurve(float [] XValue,float [] YValue,Pen p)
{

int i;
float Y;
this.max=YValue[0];
for(i=0;i<YValue.Length;i++)
if (YValue[i]>max)
this.max=YValue[i];//求出数组YValue中的最大值

//将实际值与刻值进行换算

int ScaleX=YValue.Length;
this.ScaleCount=ScaleX; 

//处理X轴刻度
int XDate=this.x-20;//实际画出的X轴长
XDate-=20;//留下20的长度不作处理

int SX=XDate/(ScaleX);//求出X轴单个刻度的长
int StartX=30;//开始时的X轴长

//处理Y轴的刻度
int YRen=this.y-20;//实际Y轴长度
YRen-=20;//留下20的长度不作处理
int SY=YRen/15;//Y轴固定画15个刻度，每个刻实际长度
int StartY=this.y-10;//开始时Y轴的开始点;if (max<150)

if (this.max<150)
{
int EndX=0;
int EndY=0;
for(i=30;i<this.x+10;i++)
{
Y=MakeEquation((float)i/SX,XValue,YValue);

int CurveX=i;

//计Y点的开始位置
float ScaleYValue=Y/10;//求对应Y轴的标筌刻度

float ScaleY=ScaleYValue*SY;//求对应Y轴的实际长度
int CurveY=StartY-(int)ScaleY;//得出Y轴的实际刻度
//int CurveY=(int)Y;

if (i==30)//第一次执行时不画点，只对其进行赋开始点位置
{
EndX=CurveX;
EndY=CurveY;
continue;
}
Point P1=new Point(EndX,EndY);//定义开始点
Point P2=new Point(CurveX,CurveY);//定义结束点
this.g.DrawLine(p,P1,P2);

EndX=CurveX;//将开始的点作为下次的起始点
EndY=CurveY;

}
}
}

//以下部份为 利用正交多项式求最小二乘解 算法

/// <summary>
/// 求最小二乘法
/// </summary>
/// <param name="X"></param>
/// <param name="x"></param>
/// <param name="y"></param>
/// <returns></returns>
private float MakeEquation(float X,float [] x,float [] y)
{
int i,j;
float a0;
float a1,a11;
float a2,a21,b1;
float Y;//返回值
j=x.Length;//得到数据的个数，即曲线拟合公式中的m

float [] y0=new float[j];
float [] y1=new float[j];
float [] y2=new float[j];

//求出系数a0的值 对应值y0恒等于1
for(i=0;i<j;i++)
y0[i]=1;//得到y0的值
a0=ExpA(y0,y);//得到a0的值

//求出系数a1的值和对应值y1
a11=ExpA(y0,x);//得到a11的值
for(i=0;i<j;i++)
y1[i]=x[i]-a11;//得到y1的值
a1=ExpA(y1,y);//得到a1的值

//求出系统a2和对应值y2
a21=ExpEXY2(x,y1)/ExpEX2(y1);
b1=ExpEX2(y1)/ExpEX2(y0);//得到b1的值
for(i=0;i<j;i++)
y2[i]=(x[i]-a21)*(x[i]-a11)-b1;//得到y2的值
a2=ExpA(y2,y);

Y=a0*1+a1*(X-a11)+a2*((X-a21)*(X-a11)-b1);//据公式y=a0*y0(X)+a1*y1(X)+a2*y2(X)
return Y;
}

/// <summary>
/// 求平方函数 公式： x*x
/// </summary>
/// <param name="x"></param>
/// <returns></returns>
private float Exp2(float x)
{
return x*x;
}
/// <summary>
/// 求公式A的值，即求多项式的系数ai的值 公式： Ex[i]*y[i]/Ex[i]*x[i]
/// </summary>
/// <param name="x"></param>
/// <param name="y"></param>
/// <returns></returns>
private float ExpA(float [] x,float [] y)
{
return ExpEXY(x,y)/ExpEX2(x);
}

private float ExpEXY2(float [] x,float [] y)
{
int i;
float sum=0;
for(i=0;i<x.Length;i++)
sum+=x[i]*Exp2(y[i]);
return sum;
}
/// <summary>
/// 求平方后的连加值 公式：Ex[i]*x[i]
/// </summary>
/// <param name="x"></param>
/// <returns></returns>
private float ExpEX2(float [] x)
{
int i;
float sum=0; 
for(i=0;i<x.Length;i++)
{
sum+=Exp2(x[i]);
}
return sum;
}

/// <summary>
/// 求X*Y后的连加值 公式：Ex[i]*y[i]
/// </summary>
/// <param name="x"></param>
/// <param name="y"></param>
/// <returns></returns>
private float ExpEXY(float [] x,float [] y)
{
int i;
float sum=0;
for (i=0;i<x.Length;i++)
{
sum+=x[i]*y[i];
}
return sum;
}
}
}