[C言語]最小二乗法を用いて数列を一次式で近似するプログラムを作ってみた[コード付き]

今回は最小二乗法を扱ってみます
なお今回の記事はコードの実装のみで、
実際に動かしてみるのは次回の記事です。

(追記　次回の記事です)

今回やりたいことを説明
最小二乗法とは
コードの実装
計算の部分は完成

今回やりたいことを説明

何をやりたいかって言うと

こんな感じでグラフ上にいくつか点があって、
これらの点すべての関係を表す関数を描けと言われたとします。

つまり尤もらしい関数を描けってことね

上と下の関数だと間違いなく下の関数の方がグラフ上の点と相関がありますよね。
今回はグラフ上に点が与えられた時にその関数を最小二乗法を用いて見つけるプログラムを組みます。

最小二乗法とは

ここは数学のブログではないので最小二乗法自体の証明などは省きます。
使えればヨシって感じ。

まずN個の数字を持つ数列bをy=f(z)に近似する場合を考えます
もう少し具体的に説明すると

このように近似するときは係数のx1とx2を求めるということです。

いま、ｚは連続の値なので
z=ntとして離散化します。(そうしないとプログラムで扱えない)
そこから

このようにJを求めます。
このJが最小になるときのx1とx2が求めるものです。これを最小二乗法と呼びます。

おお.. 原理はわからんけどすげぇ
でもどうやってそのJを求めるんだろう

ちゃんと方法があります。
求めるものが係数(x1,x2)なら

この等式を正規方程式と呼びます。
これを解けば係数が分かるんですね、不思議。

コードの実装

つまり今回は上の正規方程式が解ければ良いってことです。
なのでこの方程式を解く関数を作っていきます。

まず左辺のAtAの計算、
それから右辺のAtbを計算します。

この時点でA’x=b’の形になってますね。
ここの部分の計算は前回のガウスの消去法を用いた掃き出し法で計算します。

このプログラムでは、
係数行列と定数項ベクトルが入ったファイルを渡せば、解を計算してくれます。
（この記事も読んでくれると嬉しいです）

つまり今回の場合は、
A'(AtA)が係数行列
b'(Atb)が定数項ベクトルです。

なので前回のプログラムに渡せるようにファイルに格納する関数も作ります。

正規方程式の行列Aを作成する関数

//正規方程式で用いる行列Aの作成　1次式
void CREA1(int N, double A[N][2]){
    double t = 1/(double)N;
    for(int i=0; i<N; i++){
        A[i][0]=1;
        A[i][1]=t*(i+1);
    }
}

仮引数のNは数列の要素の数、配列Aはそのまま行列Aですね
行列Aの通りにforループを使って配列に値を入れていきます。
またtは１/Nとします。
Nがint型なのでしっかり(double)でキャストしておきました。

正規方程式の左辺を計算する関数

//正規方程式の左辺の計算
void CALA(int N, int c, double A[N][c], double L[c][c]){
    //転置したAと普通のAの積
    for(int n=0; n<c; n++){
        for(int m=0; m<c; m++){
            double S=0;
            for(int k=0; k<N; k++){
                S = S + A[k][n]*A[k][m];
            }
                L[n][m] = S; //Lに結果を入れる
        }　
    }
}

仮引数のcは係数の数です。
後に２次関数にも使えるようにするためにそうしました。
とりあえず今回は１次式なので、渡すときはc=2とします。

Lは計算結果を格納する配列です。

forループではAとAの転置の積を添字に気をつけながら計算してます。

正規方程式の右辺を計算する関数

void CALB(int N, int c, double A[N][c], double B[N], double M[c]){
    for(int i=0; i<c; i++){
        double S=0;
        for(int k=0; k<N; k++){
            S = S + A[k][i]*B[k];
        }
        M[i] = S; //Mに結果を入れる
    }
}

同じように仮引数のcは係数の数です。

こちらも同じようにforループで積を計算して結果を配列Mに格納します。

行列式をファイルに格納する関数

void INP(int c, double L[c][c], double M[c]){
    int i,j;
    FILE *fp;
    
    if((fp = fopen("ans1.txt","w"))==NULL)
        printf("ファイルを開けません\n");
    else{
        for(i=0; i<c; i++){
            for(j=0; j<c; j++){
                fprintf(fp,"%f ",L[i][j]);
            }
            fprintf(fp,"%f\n",M[i]);
        }
        fclose(fp);
    }
}

左に係数行列、
右に定数項ベクトルの形でファイルに書き出します。

あとは以下のコードで行列のサイズを2とし今作ったファイルを渡せば答えの係数が返ってきます。

void SHOW_2(int N, double A[N][N]){
    int i, j;
    for(i=0; i<N; i++){
        for(j=0; j<N; j++){
            printf("%f ",A[i][j]);
        }
        printf("\n");
    }
}

//1次元配列の中身を表示する
void SHOW_1(int N, double B[N]){
    int i;
    for(i=0; i<N; i++){
        printf("%f\n",B[i]);
    }
}

//ピボットの列の中で最も大きい値を持つ場所の添字を返す
int PIBOT(int N, double A[N][N], int k){
    double D[N];
    int p=0;
    int l;
    
    for(l=0; l<N; l++)
        D[l] = 0;
    
    
    for(l=k; l<N; l++){
        D[l] = A[l][k]*A[l][k];//絶対値で最大値を探すので二乗する
    }
    
    //最大値を探して保存
    double m = D[0];
    for(int a=1; a<N; a++)
        if(m<D[a]){
            m=D[a];
        }
    
    if(m==0){
          printf("解が一意に定まりません");
          exit(0);//exit関数を用いてきプログラムを終了する
      }
      
    
    for(int b=0; b<N; b++)
        if(m==D[b])
            p = b;//保存した最大値と同じ値をもつ添字を保存する
    
    
    return p;
}

//k列において、k行とPIBOTで得た行の交換
void CHANGE(int N, double A[N][N], double B[N], int k, int p){
    int i;
    double S;
    double T;
    
    T = B[k];
    B[k] = B[p];
    B[p] = T;
    
    for(i=0; i<N; i++){
        S = A[k][i];
        A[k][i] = A[p][i];
        A[p][i] = S;
    }
    
}

void DASH(int N, double A[N][N], double B[N], int k, int i){
    int l;
    double r = A[i][k]/A[k][k];
    
    for(l=k; l<N; l++){
        A[i][l] = A[i][l] - r*A[k][l];
    }
    
    B[i] = B[i] - r*B[k];
}


//ガウスの消去法による三角化
void DELETE(int N, double A[N][N], double B[N]){
    int k;
    int p=0;
    for(k=0; k<N-1; k++){
        p = PIBOT(N,A,k);
        CHANGE(N, A, B, k, p);
        for(int i=k+1; i<N; i++){
            DASH(N, A, B, k, i);
        }
    }
}

//ファイルから行列AとBを取り出す
void SCAN(int N, double A[N][N], double B[N], char s[NMAX]){
    double b = 0;
    FILE *fp;
    
    
    if((fp = fopen(s,"r"))==NULL)
        printf("ファイルを開けません\n");
    else{
        for(int i=0; i<N; i++){
            double a[N];
            for(int j=0; j<N; j++){
                fscanf(fp,"%lf",&a[j]);
                A[i][j] = a[j];
            }
            fscanf(fp,"%lf",&b);
                  B[i] = b;
        }
        fclose(fp);
    }
}

//解を公式から求めて配列に格納
void ANS(int N, double A[N][N], double B[N], double X[N]){
    int i,j;
    X[N-1] = B[N-1]/A[N-1][N-1];
    
    for(i=N-2; i>=0; i--){
        
        double S=0;
        for(j=i+1; j<N; j++)
            S = S + A[i][j] * X[j];
        
        X[i] = (B[i]-S)/A[i][i];
    }
    
}

int main(void) {
    int N;
    printf("行列のサイズを教えて下さい:");
    scanf("%d",&N);
     
    double A[N][N];
    double B[N];
    double X[N];
     
    
    char s[NMAX];
    printf("使用するファイルを入力してください:");
    scanf("%s",s);
    
    Q(N,s);
    SCAN(N,A,B,s);
     
    SHOW_2(N,A);
    printf("\n");
    SHOW_1(N,B);
    printf("\n");
    printf("三角化\n");
    DELETE(N,A,B);
     
    SHOW_2(N,A);
    printf("\n");
    SHOW_1(N,B);
    printf("\n");
    
    ANS(N,A,B,X);
    printf("解\n");
    SHOW_1(N,X);
  
    return 0;
}

これが上でも紹介したガウスの消去法による連立方程式を解くプログラムです。
詳しくは以下の記事で解説しているのでよければ確認してくださいorz