import java.util.Scanner;
public class Least_square_fit {
public static double Least_square_method(int n,int m,double X[],double Y[],double A[],double err[],double sum[],double my_sum,double bel[],double alp[]){
double S1[]=new double[m+1];//S1存放前一次多项式的值，范围为S1[0]~S1[m]
double S0[]=new double[m+1];//S0存放前两次多项式的值，范围为S0[0]~S0[m]
double SS[]=new double[m+1];//用于交换
double AU=0,AL=0,alp_L=0,alp_U=0,bel_U=0,bel_L=0;//AU为计算A的分子，AL为计算A的分母，alp_L、alp_U分别为计算alpha的分母和分子，alp_L_bel为计算belta的分母
double sum_temp[]=new double[m+1];///////////
double error=0;//误差的平方和
my_sum=0;
double my_sumtemp=0;
boolean flag=true;
/*计算A[0],alp[1]*/
for(int i=0;i<=m;i++)
{
AU+=Y[i];
AL++;
alp_L++;
alp_U+=X[i];
S0[i]=1;
}
A[0]=AU/AL;
bel_L=AL;
alp[1]=alp_U/alp_L;
my_sum+=A[0]*1;
for(int i=0;i<=m;i++){
sum[i]+=A[0]*1;//////////////////
}
/*计算A[1],alp[2],bel[1]*/
AU=0;AL=0;alp_L=0;alp_U=0;bel_U=0;//变量清零
double temp=0;
for(int i=0;i<=m;i++)
{
temp=(X[i]-alp[1]);
S1[i]=temp;
AU+=Y[i]*temp;
AL+=temp*temp;
alp_U+=X[i]*temp*temp;
}
alp_L=AL;
A[1]=AU/AL;
alp[2]=alp_U/alp_L;
bel_U=AL;
bel[1]=bel_U/bel_L;
my_sum+=A[1]*(X[1]-alp[1]);
for(int i=0;i<=m;i++){
sum[i]+=A[1]*(X[i]-alp[1]);///////////
}
/*递推计算A[2]~A[n-1],alp[3]~alp[n],bel[2]~bel[n-1]*/
for(int j=3;j<=n;j++){
AU=0;AL=0;alp_L=0;alp_U=0;bel_U=0;bel_L=0;//每次计算变量清零
for(int ii=0;ii<=m;ii++){
SS[ii]=S1[ii];
}
for(int i=0;i<=m;i++){
sum_temp[i]=(X[i]-alp[j-1])*S1[0]-bel[j-2]*S0[0];///////////////////
}
for(int i=0;i<=m;i++){
if(flag){
my_sumtemp=(X[1]-alp[j-1])*S1[0]-bel[j-2]*S0[0];
}
bel_L=bel_L+S1[i]*S1[i];
S1[i]=(X[i]-alp[j-1])*S1[i]-bel[j-2]*S0[i];
alp_L=alp_L+S1[i]*S1[i];
alp_U=alp_U+X[i]*S1[i]*S1[i];
S0[i]=SS[i];
AU=AU+Y[i]*S1[i];
flag=false;
}
flag=true;
bel_U=alp_L;
AL=alp_L;
alp[j]=alp_U/alp_L;
bel[j-1]=bel_U/bel_L;
A[j-1]=AU/AL;
my_sum+=A[j-1]*my_sumtemp;
for(int i=0;i<=m;i++)
{
sum[i]+=A[j-1]*sum_temp[i];
}
}
/*计算A[n]*/
AU=0;AL=0;alp_L=0;alp_U=0;bel_U=0;bel_L=0;//变量清零
for(int ii=0;ii<=m;ii++){
SS[ii]=S1[ii];
}
flag=true;
for(int i=0;i<=m;i++){
sum_temp[i]=(X[i]-alp[n])*S1[0]-bel[n-1]*S0[0];
}
for(int i=0;i<=m;i++){
if(flag){
my_sumtemp=(X[1]-alp[n])*S1[0]-bel[n-1]*S0[0];
}
S1[i]=(X[i]-alp[n])*S1[i]-bel[n-1]*S0[i];
AL=AL+S1[i]*S1[i];
S0[i]=SS[i];
AU=AU+Y[i]*S1[i];
flag=false;
}
A[n]=AU/AL;
my_sum+=A[n]*my_sumtemp;
for(int i=0;i<=m;i++){
sum[i]+=A[n]*sum_temp[i];
}
/*返回误差的平方和*/
for(int i=0;i<=m;i++)
{
err[i]=sum[i]-Y[i];
error+=err[i]*err[i];
}
return error;
}
public static void main(String[] args) {
Scanner scan=new Scanner(System.in);
System.out.println("输入多项式次数：");
int n=scan.nextInt();
System.out.println("输入数据点个数：");
int mm=scan.nextInt();
int m=mm-1;//为方便观察使用m
double bel[]=new double[n];//系数belta从bel[1]~bel[n-1]
double alp[]=new double[n+1];//系数alpha从alpha[1]~alpha[n]
double X[]=new double[m+1];//X存放mm个数据点横坐标值，范围为X[0]~X[m]
double Y[]=new double[m+1];//Y存放mm个数据点纵坐标值，范围为Y[0]~Y[m]
double A[]=new double[n+1];//A存放多项式系数，范围A[0]~A[n]
double error;//误差平方和
double err[]=new double[m+1];//记录个点误差
double sum[]=new double[m+1];//记录个点拟合值
double my_sum=0;
System.out.println("输入X:");
for(int i=0;i<=m;i++){
X[i]=scan.nextDouble();
}
System.out.println("输入Y:");
for(int i=0;i<=m;i++){
Y[i]=scan.nextDouble();
}
error=Least_square_method(n, m, X, Y,A,err,sum,my_sum,bel,alp);
System.out.println("多项式系数分别为：");//输出多项式系数
for(int i=0;i<=n;i++){
System.out.print("A["+i+"]="+A[i]+" ");
}
System.out.println();
System.out.println("alpha为：");//输出系数alpha
for(int i=1;i<=n;i++){
System.out.print("alp["+i+"]="+alp[i]+" ");
}
System.out.println();
System.out.println("belta为：");//输出系数belta
for(int i=1;i<=n-1;i++){
System.out.print("bel["+i+"]="+bel[i]+" ");
}
System.out.println();
/*输出误差相关*/
System.out.println("各点拟合值为：");
for(int i=0;i<=m;i++)
{
System.out.println(sum[i]+" ");
}
System.out.println();
System.out.println("各点误差为：");
for(int i=1;i<=m;i++){
System.out.print(err[i]+" ");
}
System.out.println();
System.out.println("误差平方和为：");
System.out.println(error);
System.out.println("my_Sum:"+my_sum);
}
}
http://www.oschina.net/code/snippet_574827_43214
