[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);
    }
}

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

計算の部分は完成

とりあえず今回のコードと、
前回のガウスの消去法のコードがあれば一次式に近似できますね

次回は実際に問題を作ってみて近似できているか検証してみます。

ではまた〜

(追記　次回の記事です)