[C言語]最小二乗法を用いて数列を一次関数に近似してgnuplotで表示してみた[コード付き]

今回はgnuplotを用いて、前回の最小二乗法のプログラムが上手くいけてるか検証してみます。


前回の記事はこちら

また今回は全体のコードを下に貼っておくので参考にどうぞ。



検証方法

まずは問題となる数列を作りましょう。
数列というよりは二次元座標上の点ですね。

完全ランダムな点をNを発生させてその近似の関数を求めてもよいのですが、
完全にバラバラの点だとうまく近似できているか判断できません。



なので今回は答えとなる関数を決めておきます。


その関数をxが0~1の範囲内でN等分して、
N個の点に離散化します。

そしてそのN個の点に雑音を付加します。

雑音というのは簡単に言えばランダムな値ですね。
そうすると上の図のオレンジの点が

この図の青色の点のようになります。

この青色の点に前回のプログラムをかけて
最初に決めた答えとなる関数に近似できるか検証します。



数列を作る関数を実装

N個の点に付加する雑音ですが、
以前に作ったガウス性雑音を用います。

ガウス性雑音自体は以前に紹介したので気になる方はどうぞ

コードです

//数列を生成しファイルに入力する
//ave,verは雑音における平均と分散
void CRE1(int N, double x1, double x2, double ave, double ver){
    FILE *fp;
      
    if((fp=fopen("kk.txt","w"))==NULL){
        printf("ファイルに書き込めません");
    }else{
        for(int i=1; i<=N; i++){
              
            double Y=0;
            Y=x1+x2*((double)i/N)+ave+ver*BOX();//1次式にガウス性雑音を付加
          
            fprintf(fp,"%f\n",Y);
              
        }
        fclose(fp);
    }
}

まず関数BOXはガウス性雑音を返す関数です。
(BOXのコードは下の全体のコードに貼ってあるので確認どうぞ)

仮引数のNが点の個数、x1,x2はそれぞれ係数
aveが雑音の平均でverが分散です。

forループの中で

Y=x1+x2*((double)i/N)+ave+ver*BOX();

このように計算してます。
i/Nとすることで、iが1~Nと動くので0~1の範囲で離散化できてます。
そこにBOXでガウス性雑音を足してます。

またファイルはここでは”kk.txt”としてますが別に何でも良いです。



検証

数列のみをgnuplotでプロット

y = a + bxにおいて
a=1、b=2
発生する点が100個で離散化
ガウス性雑音の平均0、分散0.5で数列を作ってみました

うっすーらy=1+2xって感じがしますね。



行列式を求める

この数列を用いて、前回のプログラムを掛けて行列式を求めると

こんな行列式が返ってきました。



ガウスの消去法でこの行列式を解く

以前に作ったガウスの消去法のプログラミングに上の行列式を渡して解かせてみると、



解が
1.017196
1.927613
ということは
1
2
にまあまあ近似できてますね

gnuplotで点と一緒にグラフに描くと


いい感じじゃないですか。






雑音と分散を変えた場合の結果

雑音と平均以外はそのままで何個か検証してみました。

平均0 分散0.3

平均0.3 分散0.3

平均0 分散0(雑音を与えて無いのと等しい)


大丈夫そうですね。
どのパターンでもしっかり近似できてます。



全体のコード

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
//ボックスミュラー法によるガウス性雑音の生成
double BOX(void){
    double y1,y2,x1;
    y1 = (double)rand()/(double)RAND_MAX;
    y2 = (double)rand()/(double)RAND_MAX;
              
    x1 = sqrt(-2*log(y1))*cos(2*M_PI*y2);
    
    return x1;
}

//数列を生成しファイルに入力する
//ave,verは雑音における平均と分散
void CRE1(int N, double x1, double x2, double ave, double ver){
    FILE *fp;
      
    if((fp=fopen("kk.txt","w"))==NULL){
        printf("ファイルに書き込めません");
    }else{
        for(int i=1; i<=N; i++){
              
            double Y=0;
            Y=x1+x2*((double)i/N)+ave+ver*BOX();//1次式にガウス性雑音を付加
          
            fprintf(fp,"%f\n",Y);
              
        }
        fclose(fp);
    }
    
}


//正規方程式で用いる行列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);
    }
}


//数列が入ったファイルからプログラム内の配列に出力する
void SCAN1(int N, double B[N]){
    FILE *fp;
    
    if((fp=fopen("kk.txt","r"))==NULL){
        printf("ファイルを開けません");
    }else{
        for(int i=0; i<N; i++){
            double b=0;
            fscanf(fp,"%lf",&b);
            B[i] = b;
        }
        fclose(fp);
    }
}

//正規方程式の左辺の計算
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に結果を入れる
          
        }
    }
}

//正規方程式の右辺の計算
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に結果を入れる
    }
}

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

int main(void) {
   int N;
    int x1,x2;
    printf("発生する点の数を入力してください:");
    scanf("%d",&N);
    printf("近似する一次式の係数を入力してください\nx1:");
    scanf("%d",&x1);
    printf("\nx2:");
    scanf("%d",&x2);
    
    double ave,ver;
    printf("発生する雑音の平均を入力してください:");
    scanf("%lf",&ave);
    printf("発生する雑音の分散を入力してください:");
    scanf("%lf",&ver);
    
    double A[N][2];
    double B[N];
    
    CRE1(N,x1,x2,ave,ver);
    SCAN1(N,B);
    CREA1(N,A);
    
    double L[2][2];
    double M[2];
       
    CALA(N,2,A,L);
    CALB(N,2,A,B,M);
       
    INP(2,L,M);
    
    return 0;
}



このコードで行列式の作成までが完了です。

この記事ではそこからは、以前に紹介したガウスの消去法を用いた掃き出し法のプログラミングの一部を利用してます。
よかったらそちらも参考にどうぞ

ガウスの消去法後編

今回はここまでです。


ではまた〜。

syumokuzame

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

コメントする