采用最小二乘法拟合曲线：请大家帮忙。(100分)

我想通过一些坐标点采用最小二乘法来拟合曲线，次数最好大于3次。我找了一些程序，可是不知道怎样调用函数（例如下面的），哪位大侠能举个例子告诉我怎样用吗？？
***************************************************************
* 本算法用最小二乘法依据指定的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);
}

我用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>max)
this.max=YValue;//求出数组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/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 = 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>max)
this.max=YValue;//求出数组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=1;//得到y0的值
a0=ExpA(y0,y);//得到a0的值

//求出系数a1的值和对应值y1
a11=ExpA(y0,x);//得到a11的值
for(i=0;i<j;i++)
y1=x-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=(x-a21)*(x-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*y/Ex*x
/// </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*Exp2(y);
return sum;
}
/// <summary>
/// 求平方后的连加值 公式：Ex*x
/// </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);
}
return sum;
}
/// <summary>
/// 求X*Y后的连加值 公式：Ex*y
/// </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*y;
}
return sum;
}
}
}

赚钱的好方法，绝对真实！！！
事先声明，以下以讨论为主，不喜欢的不必参与
　　现在大部分的网络兼职，无非就是帮网络广告商挂广告条，靠积分赚钱。而大部分人都对此不屑一顾，其中也包括我自己。一个月前，在极度无聊的情况下，我尝试了两家国内较为知名的网络广告商：NewsBar和八趣通宝。这两家的机制差不多，操作也差不多，无非就是上网的时候打开广告条，然后最小化让它自己记点就好了。上个月17号和21号，我分别收到了NewsBar和八趣通宝116.3＋121.5共计237.8元的银行转帐。因为是刚开始做，钱不多，主要靠下线发展，就我知道的，有人单单做NewsBar，一个月最高做了1700多，我是不敢奢望了，以后一个月有四五百就满足了，毕竟是免费且无需多少精力投入的东东，赚多少算多少吧。
　　
绝对不要相信那些网赚可以月入万元的广告，纯属放屁，真正的网赚，你做单个的话，一个月有七八百已 经很多了，当然可以几个一起做，但一定要找可靠的广告商，国外的就免谈了，汇款太繁琐，还是国内直接银行转帐，真假立辨。但是，刚开始做的时候，收入是很少的，等到下线足够多的时候，自己只要做足每月2000点（五六天时间）就可以有一笔不小的收入了，重要的是要有耐心，一开始确实是比较少的：）而且作为下线你也不用担心，给上线的奖励是额外的，并不会侵犯你本人的劳动成果
　　
　　提供三个连接，有兴趣的朋友可以试试，没兴趣的也不必骂人，飘走就是了
　　
　　1、NewsBar，知名广告商，运作稳定，推荐
　　
www.ads4cn.com/newsbar/refferer.asp?jay1717　
　　2、八趣通宝，知名广告商，运作稳定，推荐
www.8qu.net/register.asp?net=jay1717　　
　　
　　3、soho广告，完全模仿NewsBar，目前系统尚不稳定，但积分很快，
www.sohoads.com/sabar/reg.asp?sj=jay1717
当然，还是要鼓励大家一起挂，毕竟自己挂的也是钱啊！ 从今天开始，记住了：如果开QQ, 玩棋牌游戏, 逛论坛的时候也开这三个小软件,一边娱乐一边赚钱,一举两得~~~~~~~~~~简单赚钱