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