用正交多项式作最小二乘拟合的java实现（转）

柔情只为你懂 2022-03-30 05:25 177阅读 0赞

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