C言語で多項式回帰アルゴリズムを計算する方法【サンプルコード付き】

回帰分析とは、従属変数と独立変数の間の関係を明らかにするための予測モデリング手法の一つです。過去のデータに基づいて将来の値を予測したり、変数間の影響の強さを評価したりする際に、統計解析や機械学習の分野で広く活用されています。

多項式回帰とは

多項式回帰は、回帰分析の一種であり、独立変数 x と従属変数 y の関係を、x の n 次多項式としてモデル化する手法です。単純な直線（一次回帰）では表現しきれない、曲線的な複雑な傾向を捉えることができるため、実務のデータ分析でも多くの場面で利用されています。

プログラムの仕組み

以下のCプログラムは、最小二乗法に基づいて多項式回帰の係数を求めます。処理の流れは次のとおりです。

1. 多項式の次数 m とデータ数 n を入力する
2. 各データ点 (x, y) の値を入力する
3. 観測データから正規方程式の係数行列と右辺ベクトルを構築する
4. ガウスの消去法を用いて連立方程式を解く
5. 多項式の各係数 a[1] 〜 a[m+1] を出力する

サンプルコード

以下は、多項式回帰アルゴリズムを計算するCプログラムです。

#include<math.h>
#include<stdio.h>
#include<conio.h>
main(){
   int i,j,k,m,n;
   float x[20],y[20],u,a[10],c[20][20],power,r;
   printf("enter m,n:");
   scanf("%d%d",&m,&n);
   for(i=1;i<=n;i++){
      printf("enter values of x and y");
      scanf("%f%f",&x[i],&y[i]);
   }
   for(j=1;j<=m+1;j++)
   for(k=1;k<=m+1;k++){
      c[j][k]=0;
      for(i=1;i<=n;i++){
         power=pow(x[i],j+k-2);
         c[j][k]=c[j][k]+power;
      }
   }
   for(j=1;j<=m+1;j++){
      c[j][m+2]=0;
      for(i=1;i<=n;i++){
         r=pow(x[i],j-1);
         c[j][m+2]=c[j][m+2]+y[i]*r;
      }
   }
   for(i=1;i<=m+1;i++){
      for(j=1;j<=m+2;j++){
         printf("%.2f\t",c[i][j]);
      }
      printf("\n");
   }
   for(k=1;k<=m+1;k++)
      for(i=1;i<=m+1;i++){
         if(i!=k){
            u=c[i][k]/c[k][k];
            for(j=k;j<=m+2;j++){
               c[i][j]=c[i][j]-u*c[k][j];
         }
      }
   }
   for(i=1;i<=m+1;i++){
      a[i]=c[i][m+2]/c[i][i];
      printf("a[%d]=%f\n",i,a[i]);
   }
   getch();
}

実行結果

上記のプログラムをコンパイルして実行すると、次のような結果が出力されます。

enter m,n:4 5
enter values of x and y1 1
enter values of x and y2 3
enter values of x and y1 2
enter values of x and y1 2
enter values of x and y1 1
5.00  6.00  8.00  12.00  20.00  9.00
6.00  8.00  12.00 20.00  36.00  12.00
8.00  12.00 20.00 36.00  68.00  18.00
12.00 20.00 36.00 68.00  132.00 30.00
20.00 36.00 68.00 132.00 260.00 54.00
a[1]=1.750000
a[2]=-2.375000
a[3]=2.000000
a[4]=0.500000
a[5]=-0.375000

結果の見方

出力の前半に表示されている行列は、観測データから構築された正規方程式の拡大係数行列です。そして最後に表示される a[1] 〜 a[5] の値が、求める多項式の各係数に対応します。

この例では、4次多項式 y = a₁ + a₂x + a₃x² + a₄x³ + a₅x⁴ に対して、y = 1.75 − 2.375x + 2x² + 0.5x³ − 0.375x⁴ という回帰式が得られます。この式に x の値を代入することで、対応する y の予測値を計算できます。