最小２乗法により、ｎ次多項式で近似を行うプログラムである。今回は３次式で近似を行ってみる。プログラム中で近似を行ったデータは、適当に決めたものである。

またガウスの消去法のプログラムはガウスの消去法の数値計算例(Java)のプログラムを少し変えたものである。

import java.io.*;

class Co
{
        static final int N=4;   //ガウスの消去法における未知数の数（下のnと同じ値にすること）
        static final int n=4;   //n-1次の多項式で近似（１以上）
        static final int S=5;   //データの個数
}

class Saisyounizyouhou
{
        //配列xxは解、すなわち多項式の係数を入れるためのもの
        static double xx[]=new double[Co.N];
        static double a[][]=new double[Co.n][Co.n+1];

        public static void main(String arg[])
        {
                double x[]={1.0,2.0,3.0,4.0,5.0},y[]={7.987,2.986,1.998,2.224,5.678};

                Saisyounizyouhou sn=new Saisyounizyouhou();

                /*データの横軸に値する数値を配列xに小さい順に入れ、
                 それぞれの横軸の値における縦軸の値を配列yに入れ関数saiに渡す*/
                sn.sai(x,y);
        }

        void sai(double x[],double y[])
        {
                Saisyounizyouhou sn=new Saisyounizyouhou();

                int i,j,k;
                double X,Y;

                try
                {
                        FileWriter fw1=new FileWriter("output1.data");
                        FileWriter fw2=new FileWriter("output2.data");

                        /*配列aの初期化として0を代入しておく*/
                        for(i=0;i<Co.n;i++) 
                        {
                                for(j=0;j<Co.n+1;j++) 
                                {
                                        a[i][j]=0.0;
                                }
                        }

                        /*ガウスの消去法で解く行列の作成*/
                        for(i=0;i<Co.n;i++) 
                        {
                                for(j=0;j<Co.n;j++) 
                                {
                                        for(k=0;k<Co.S;k++) 
                                        {
                                                a[i][j]+=Math.pow(x[k],i+j);
                                        }
                                }
                        }
                        for(i=0;i<Co.n;i++) 
                        {
                                for(k=0;k<Co.S;k++) 
                                {
                                        a[i][Co.n]+=Math.pow(x[k],i)*y[k];
                                }
                        }
                        /*ガウスの消去法の実行*/
                        sn.gauss();

                        /*GNUPLOTで表示するために最小２乗法による関数のデータをファイル保存*/
                        for(X=x[0]-10.0;X<x[Co.S-1]+10.0;X+=0.01) 
                        {
                                Y=0.0;
                                for(i=0;i<Co.N;i++) 
                                {
                                        Y+=xx[i]*Math.pow(X,i);
                                }
                                fw1.write(X + " " + Y + "\n");
                        }

                        /*GNUPLOTで表示するために、最小２乗法に使われたデータを保存*/
                        for(i=0;i<Co.S;i++) 
                        {
                                fw2.write(x[i] + " " + y[i] + "\n");
                        }

                        fw1.close();
                        fw2.close();
                }
                catch(Exception e)
                {
                        System.out.println(e);
                }
        }

        void gauss()
        {
                double x[]=new double[Co.N];
                double b[][]=new double[1][Co.N+1];

                int i,j,k,l,pivot;
                double p,q,m;

                for(i=0;i<Co.N;i++) 
                {
                        m=0;
                        pivot=i;

                        for(l=i;l<Co.N;l++) 
                        {
                                if(Math.abs(a[l][i])>m) //i列の中で一番値が大きい行を選ぶ 
                                {   
                                        m=Math.abs(a[l][i]);
                                        pivot=l;
                                }
                        }

                        if(pivot!=i)    //pivotがiと違えば、行の入れ替え 
                        {                          
                                for(j=0;j<Co.N+1;j++) 
                                {
                                        b[0][j]=a[i][j];        
                                        a[i][j]=a[pivot][j];
                                        a[pivot][j]=b[0][j];
                                }
                        }
                }

                for(k=0;k<Co.N;k++) 
                {
                        p=a[k][k];              //対角要素を保存
                        a[k][k]=1;              //対角要素は１になることがわかっているから

                        for(j=k+1;j<Co.N+1;j++) 
                        {
                                a[k][j]/=p;
                        }

                        for(i=k+1;i<Co.N;i++) 
                        {
                                q=a[i][k];

                                for(j=k+1;j<Co.N+1;j++) 
                                {
                                        a[i][j]-=q*a[k][j];
                                }
                                a[i][k]=0;              //0となることがわかっているところ
                        }
                }

                /*解の計算*/
                for(i=Co.N-1;i>=0;i--) 
                {
                        x[i]=a[i][Co.N];
                        for(j=Co.N-1;j>i;j--) 
                        {
                                x[i]-=a[i][j]*x[j];
                        }
                }

                /*行列が最後どうなったかみる*/
                for(i=0;i<Co.N;i++) 
                {
                        for(j=0;j<Co.N+1;j++) 
                        {
                                System.out.print(a[i][j] + " ");
                        }
                        System.out.println("");              
                }

                /*解の確認表示*/
                System.out.println("解は");
                for(i=0;i<Co.N;i++) 
                {
                        System.out.println(x[i]);
                        xx[i]=x[i];     //解の値をメソッドsaiで使うために配列xxに代入
                }
        }
}

緑色の点がデータ、赤色の線が最小２乗法で得られた直線である。