2010.11. 7.
石立 喬

Visual C++ 2010 Express の易しい使い方(18)

―――― 最小二乗法により、測定データに多項式を当てはめる ――――


概要
  実験結果をグラフに表示する時、測定点をプロットしてから、測定点に最も近い直線や二次曲線を描くことがある。これを自動的に行うプログラムを紹介する。この作業は「曲線の当てはめ」と呼ばれ、多元連立方程式を解く必要がある。

当てはめ
  「曲線の当てはめ(Curve Fitting)」とは、測定データとの誤差の二乗の和が最小になるような、与えられた形式の式(実験式)を求めることを言う。「易しい使い方(17)」で紹介した「補間」が測定データを絶対視し、それを必ず通る曲線を描き、中間の値を推定したのに対し、「当てはめ」は、測定データからの誤差を最小にする式を求める。
  一番単純な当てはめは直線(一次式)で近似するもので、その直線は「回帰直線」と呼ばれ、統計で良く用いられる。二次式や三次式を当てはめる方法もあるが、測定データは理論的に平方根、指数、冪などの形をとるべきものも多く、無理に高次の多項式に当てはめるのはあまり意味がない。

最小二乗法による多項式の当てはめ
 当てはめは下記の手順による。
 1) 目的とする一次式、二次式、三次式などの式を仮定する。
 2) 二乗誤差の和を最小にするために、式の各係数a、b、c、dなどについて偏微分を求め、それを0と置く。
 3) 二元(一次式の場合)、三元(二次式)、四元(三次式)などの連立方程式が得られる。
 4) 上記の連立方程式を解く。ここでは、ピボット選択付きのGauss-Jordan法を利用する。
 5) 連立方程式の解として得た各係数を用いてグラフを描く。
 6) 式と二乗平均誤差を表示する。

当てはめに用いる連立方程式
◎一次式による当てはめ(回帰直線)の場合
当てはめたい式を
  
とすると、二乗誤差の和は、データをxi、yiとし、データ点の数を n として、
  
である。これを最小にする a および b の値は、
  
で求められる。したがって、解くべき連立方程式は、
  
である。
◎二次式による当てはめの場合
 当てはめたい式を
     
とすると、途中を省略して、解くべき連立方程式は次のようになる。
  

◎三次式による当てはめの場合
 当てはめたい式を
  
とすると、解くべき連立方程式は次のようになる。
  

連立方程式を解くプログラムと使用する配列
 ここでは、一般に良く知られているGauss-Jordanの方法を用いた。なるべく精度を高め、0割りなどのエラーが生じないように、ピボット選択方式を採用した。
 三次式による当てはめ、すなわち四元連立方程式を解く場合を例として下表に示す。上の式と比較して欲しい。sx6などは、xの6乗の累計の意味である。

  a*aa[0,0] + b*aa[0,1] + c*aa[0,2] + d*aa[0,3] = aa[0,4]
  a*aa[1,0] + b*aa[1,1] + c*aa[1,2] + d*aa[1,3] = aa[1,4]
  a*aa[2,0] + b*aa[2,1] + c*aa[2,2] + d*aa[2,3] = aa[2,4]
  a*aa[3,0] + b*aa[3,1] + c*aa[3,2] + d*aa[3,3] = aa[3,4]

のように配列aa[i,j]を割り当て、aa[0,1]=sx6; ----- aa[0,4]=sx3y; ----- aa[3,3]=data_number; などと代入して、a、b、c、dを求める。
 配列のみを取り出すと、下記のような行列ができる。
 
  aa[0,0]   aa[0,1]   aa[0,2]   aa[0,3]   aa[0,4]
  aa[1,0]   aa[1,1]   aa[1,2]   aa[1,3]   aa[1,4]
  aa[2,0]   aa[2,1]   aa[2,2]   aa[2,3]   aa[2,4]
  aa[3,0]   aa[3,1]   aa[3,2]   aa[3,3]   aa[3,4]

 これを変形して

    1       0       0       0    aa[0,4]  (=b[0])
   0      1      0       0    aa[1,4]  (=b[1])
   0      0      1       0   aa[2,4]  (=b[2])
   0      0      0       1   aa[3,4]  (=b[3])

のようにするのが、Gauss-Jordan法であり、解は b[0] 〜 b[3] に得られる。
 なお、Gauss-Jordan法の原理を説明するために、プログラム中にコメント文をやや詳しく付けた。

プログラムの概要
◎「易しい使い方(17)」の補間で用いたプログラムを流用する。
 下記に示すメソッドを流用した。これらについては、説明を省略する。
  Form1_Paint()
  Form1_MouseMove()
  Form1_MouseDown()
  buttonFittingCancel_Click() ---- buttonInterCancel_Click() と名称が異なるが、内容は同一
  buttonDataCancel_Click()
◎あてはめプログラム
 1)下表の関係に従って、配列を用意し、各係数を計算する。

あてはめる式 解くべき連立方程式  各係数の配列  解を受け取る配列
一次式 二元連立方程式 aa[0,0]〜aa[1,2] bb[0]〜bb[1]
二次式 三元連立方程式 aa[0,0]〜aa[2,3] bb[0]〜bb[2]
三次式 四元連立方程式 aa[0,0]〜aa[3,4] bb[0]〜bb[3]

 2)ユーザ定義メソッドcalculateGaussJordanを呼び出し、結果をa、b、cなどに代入する。
 3)a、b、cなどを用いて、対応する直線(曲線)を描く。
 4)式を表示する。
 5)二乗誤差を計算して表示する。

プログラム

static int X0=50,X1=450;
static int Y0=60,Y1=500,Y2=540,Y3=560,Y4=580;
static int data_number=0;
static array<float>^ xdata=gcnew array<float>(20);
static array<float>^ ydata=gcnew array<float>(20);

private: System::Void buttonFirst_Click(System::Object^ sender, System::EventArgs^ e) {

   Graphics^ g=this->CreateGraphics();

   int i;
   float x1,y1,x2,y2;
   float y;
   float sx=0.0f,sy=0.0f,sx2=0.0f,sxy=0.0f;
   float error=0.0f;
   float a,b;
   array<float,2>^ aa=gcnew array<float,2>(2,3);
   array<float>^ bb=gcnew array<float>(2);

   //二元連立方程式の係数を計算
   for(i=0;i<data_number;i++){
      sx+=xdata[i];
      sy+=ydata[i];
      sx2+=xdata[i]*xdata[i];
      sxy+=xdata[i]*ydata[i];
   }
   aa[0,0]=sx2; aa[0,1]=sx;      aa[0,2]=sxy;
   aa[1,0]=sx; aa[1,1]=data_number; aa[1,2]=sy;

   //二元連立方程式を解く
   bb=calculateGaussJordan(2,aa);
   a=bb[0];
   b=bb[1];

   //直線式を描画
   x1=xdata[0];
   x2=xdata[data_number-1];
   y1=a*x1+b;
   y2=a*x2+b;
   g->DrawLine(Pens::Blue,X0+10*x1,Y0+10*(40-y1),X0+10*x2,Y0+10*(40-y2));
   if(b>=0)
      g->DrawString("一次式  y=" +a+ "*x + " +b,Font,Brushes::Blue,X0,Y2);
   else
      g->DrawString("一次式  y=" +a+ "*x - " +(-b),Font,Brushes::Blue,X0,Y2);

   //二乗誤差を表示
   for(i=0;i<data_number;i++){
      y=a*xdata[i]+b;
      error+=(float)Math::Pow(ydata[i]-y,2);
   }
   g->DrawString("二乗誤差の合計=" +error,Font,Brushes::Blue,X1,Y2);

}

private: System::Void buttonSecond_Click(System::Object^ sender, System::EventArgs^ e) {

   Graphics^ g=this->CreateGraphics();

   int i;
   float x,y;
   float xx,yy,xx_old,yy_old;
   float x1,x2,x3,x4;
   float sx=0.0f,sx2=0.0f,sx3=0.0f,sx4=0.0f,sy=0.0f,sxy=0.0f,sx2y=0.0f,error=0.0f;
   array<float,2>^ aa=gcnew array<float,2>(3,4);
   array<float>^ bb=gcnew array<float>(4);
   float a,b,c;

   //三元連立方程式の係数を計算
   for(i=0;i<data_number;i++){
      x1=xdata[i];
      x2=x1*xdata[i];
      x3=x2*xdata[i];
      x4=x3*xdata[i];
      sx+=x1;
      sx2+=x2;
      sx3+=x3;
      sx4+=x4;
      sy+=ydata[i];
      sxy+=x1*ydata[i];
      sx2y+=x2*ydata[i];
   }
   aa[0,0]=sx4; aa[0,1]=sx3; aa[0,2]=sx2;     aa[0,3]=sx2y;
   aa[1,0]=sx3; aa[1,1]=sx2; aa[1,2]=sx;      aa[1,3]=sxy;
   aa[2,0]=sx2; aa[2,1]=sx;  aa[2,2]=data_number; aa[2,3]=sy;

   //三元連立方程式を解く
   bb=calculateGaussJordan(3,aa);
   a=bb[0];
   b=bb[1];
   c=bb[2];

   //二次曲線を描画
   for(x=xdata[0];x<=xdata[data_number-1];x+=0.01f){
      y=a*x*x+b*x+c;
      xx=X0+x*10;
      yy=Y0+(40-y)*10;
      if(x==xdata[0]){
         xx_old=xx;
         yy_old=yy;
      }
      else{
         g->DrawLine(Pens::YellowGreen,xx_old,yy_old,xx,yy);
         xx_old=xx;
         yy_old=yy;
      }
   }

   //二次式を表示
   String^ string1;
   if(b>=0){
      if(c>=0) string1="二次式  y=" +a+ "x*x + " + b + "x + " + c;
      else   string1="二次式  y=" +a+ "x*x + " + b + "x - " +(-c);
   }
   else{
      if(c>=0) string1="二次式  y=" +a+ "x*x - " +(-b)+ "x + " + c;
      else   string1="二次式  y=" +a+ "x*x - " +(-b)+ "x - " +(-c);
   }
   g->DrawString(string1,Font,Brushes::Green,X0,Y3);

   //二乗誤差を表示
   for(i=0;i<data_number;i++){
      y=a*xdata[i]*xdata[i]+b*xdata[i]+c;
      error+=(float)Math::Pow(ydata[i]-y,2);
   }
   g->DrawString("二乗誤差の合計=" +error,Font,Brushes::Green,X1,Y3);

}

private: System::Void buttonThird_Click(System::Object^ sender, System::EventArgs^ e) {

   Graphics^ g=this->CreateGraphics();

   int i;
   float x,y;
   float xx,yy,xx_old,yy_old;
   float x1,x2,x3,x4,x5,x6;
   float sx=0.0f,sx2=0.0f,sx3=0.0f,sx4=0.0f,sx5=0.0f,sx6=0.0f;
   float sy=0.0f,sxy=0.0f,sx2y=0.0f,sx3y=0.0f,error=0.0f;
   array<float,2>^ aa=gcnew array<float,2>(4,5);
   array<float>^ bb=gcnew array<float>(5);
   float a,b,c,d;

   //四元連立方程式の係数を計算
   for(i=0;i<data_number;i++){
      x1=xdata[i];
      x2=x1*xdata[i];
      x3=x2*xdata[i];
      x4=x3*xdata[i];
      x5=x4*xdata[i];
      x6=x5*xdata[i];
      sx+=x1;
      sx2+=x2;
      sx3+=x3;
      sx4+=x4;
      sx5+=x5;
      sx6+=x6;
      sy+=ydata[i];
      sxy+=x1*ydata[i];
      sx2y+=x2*ydata[i];
      sx3y+=x3*ydata[i];
   }
   aa[0,0]=sx6; aa[0,1]=sx5; aa[0,2]=sx4; aa[0,3]=sx3;     aa[0,4]=sx3y;
   aa[1,0]=sx5; aa[1,1]=sx4; aa[1,2]=sx3; aa[1,3]=sx2;     aa[1,4]=sx2y;
   aa[2,0]=sx4; aa[2,1]=sx3; aa[2,2]=sx2; aa[2,3]=sx;     aa[2,4]=sxy;
   aa[3,0]=sx3; aa[3,1]=sx2; aa[3,2]=sx;  aa[3,3]=data_number; aa[3,4]=sy;

   //四元連立方程式を解く
   bb=calculateGaussJordan(4,aa);
   a=bb[0];
   b=bb[1];
   c=bb[2];
   d=bb[3];

   //三次曲線を描画
   for(x=xdata[0];x<=xdata[data_number-1];x+=0.01f){
      y=a*x*x*x+b*x*x+c*x+d;
      xx=X0+x*10;
      yy=Y0+(40-y)*10;
      if(x==xdata[0]){
         xx_old=xx;
         yy_old=yy;
      }
      else{
         g->DrawLine(Pens::OrangeRed,xx_old,yy_old,xx,yy);
         xx_old=xx;
         yy_old=yy;
      }
   }

   //三次式を表示
   String^ string1;
   if(b>=0){
      if(c>=0){
         if(d>=0) string1="三次式  y=" +a+ "x*x*x + " + b + "x*x + " + c + "x + " +d;
         else   string1="三次式  y=" +a+ "x*x*x + " + b + "x*x + " + c + "x - " +(-d);
      }
      else{
         if(d>=0) string1="三次式  y=" +a+ "x*x*x + " + b + "x*x - " +(-c)+ "x + " +d;
         else   string1="三次式  y=" +a+ "x*x*x + " + b + "x*x - " +(-c)+ "x - " +(-d);
      }
   }
   else{
      if(c>=0){
         if(d>=0) string1="三次式  y=" +a+ "x*x*x - " +(-b)+ "x*x + " + c + "x + " +d;
         else   string1="三次式  y=" +a+ "x*x*x - " +(-b)+ "x*x + " + c + "x - " +(-d);
      }
      else{
         if(d>=0) string1="三次式  y=" +a+ "x*x*x - " +(-b)+ "x*x - " +(-c)+ "x + " +d;
         else   string1="三次式  y=" +a+ "x*x*x - " +(-b)+ "x*x - " +(-c)+ "x - " +(-d);
      }
   }
   g->DrawString(string1,Font,Brushes::OrangeRed,X0,Y4);

   //二乗誤差を表示
   for(i=0;i<data_number;i++){
      y=a*xdata[i]*xdata[i]*xdata[i]+b*xdata[i]*xdata[i]+c*xdata[i]+d;
      error+=(float)Math::Pow(ydata[i]-y,2);
   }
   g->DrawString("二乗誤差の合計=" +error,Font,Brushes::OrangeRed,X1,Y4);

}

private: array<float>^ calculateGaussJordan(int n,array<float,2>^ aa){

   float p,d,coef_max,temp;
   int i,j,k;
   int pivot;
   for(k=0;k<n;k++){   //係数配列aa[i,j]の列(column)を横方向にkでスキャン
      //k列で係数の絶対値が最大となる行をさがし、pivotに設定
      coef_max=0.0;
      for(i=k;i<n;i++){   //k番目からはじめて、行(row)を縦方向にiでスキャン
         if(Math::Abs(aa[i,k])>coef_max){
            coef_max=Math::Abs(aa[i,k]);
            pivot=i;
         }
      }
      for(j=0;j<=n;j++){   //横方向にjでスキャンして、aa[k,j]とaa[pivot,j]を交換
         temp=aa[k,j];
         aa[k,j]=aa[pivot,j];
         aa[pivot,j]=temp;
      }
      //aa[k,k]が1になるように、横方向にjでスキャンして、aa[k,k]〜aa[k,n]をaa[k,k]で割る
      p=aa[k,k];
      for(j=k;j<=n;j++)
         aa[k,j]/=p;
      //列を縦方向にiでスキャンして、aa[i,k]が0になるようにする
      for(i=0;i<n;i++){
         if(i!=k){
            d=aa[i,k];
            for(j=k;j<=n;j++)
               aa[i,j]-=d*aa[k,j];
         }
      }
   }

   //解を配列bb[n]に入れて返す
   array<float>^ bb=gcnew array<float>(n);
   for(i=0;i<n;i++)
      bb[i]=aa[i,n];

   return bb;

}

プログラムの実行結果の一例
 図1において、黒丸はマウスで入力したデータ、青は「一次式」、緑は「二次式」、赤は「三次式」のボタンをそれぞれクリックした結果である。「方眼紙」の下は、各線の式と二乗誤差が表示されている。一次式から三次式に次数が高くなるほど誤差が減少している。「方眼紙」の下にはマウスが指している座標が表示されるが、この例では、マウスが「方眼紙」外にあったため「Out of Range」になっている。



図1 得られた結果の一例


「Visual C++ の勉強部屋」(目次)へ