2010.11. 7.
石立 喬
Visual C++ 2010 Express の易しい使い方（１8）
―――― 最小二乗法により、測定データに多項式を当てはめる ――――
概要
実験結果をグラフに表示する時、測定点をプロットしてから、測定点に最も近い直線や二次曲線を描くことがある。これを自動的に行うプログラムを紹介する。この作業は「曲線の当てはめ」と呼ばれ、多元連立方程式を解く必要がある。
当てはめ
「曲線の当てはめ（Curve Fitting）」とは、測定データとの誤差の二乗の和が最小になるような、与えられた形式の式（実験式）を求めることを言う。「易しい使い方（１7）」で紹介した「補間」が測定データを絶対視し、それを必ず通る曲線を描き、中間の値を推定したのに対し、「当てはめ」は、測定データからの誤差を最小にする式を求める。
一番単純な当てはめは直線（一次式）で近似するもので、その直線は「回帰直線」と呼ばれ、統計で良く用いられる。二次式や三次式を当てはめる方法もあるが、測定データは理論的に平方根、指数、冪などの形をとるべきものも多く、無理に高次の多項式に当てはめるのはあまり意味がない。
最小二乗法による多項式の当てはめ
当てはめは下記の手順による。
１） 目的とする一次式、二次式、三次式などの式を仮定する。
２） 二乗誤差の和を最小にするために、式の各係数a、b、c、dなどについて偏微分を求め、それを0と置く。
３） 二元（一次式の場合）、三元（二次式）、四元（三次式）などの連立方程式が得られる。
４） 上記の連立方程式を解く。ここでは、ピボット選択付きのGauss-Jordan法を利用する。
５） 連立方程式の解として得た各係数を用いてグラフを描く。
６） 式と二乗平均誤差を表示する。
当てはめに用いる連立方程式
◎一次式による当てはめ（回帰直線）の場合
当てはめたい式を
とすると、二乗誤差の和は、データを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法の原理を説明するために、プログラム中にコメント文をやや詳しく付けた。
プログラムの概要
◎「易しい使い方（１７）」の補間で用いたプログラムを流用する。
下記に示すメソッドを流用した。これらについては、説明を省略する。
Form1_Paint()
Form1_MouseMove()
Form1_MouseDown()
buttonFittingCancel_Click() ---- buttonInterCancel_Click() と名称が異なるが、内容は同一
buttonDataCancel_Click()
◎あてはめプログラム
１）下表の関係に従って、配列を用意し、各係数を計算する。
| あてはめる式 |
解くべき連立方程式 |
各係数の配列 |
解を受け取る配列 |
| 一次式 |
二元連立方程式 |
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] |
２）ユーザ定義メソッドcalculateGaussJordanを呼び出し、結果をa、b、cなどに代入する。
３）a、b、cなどを用いて、対応する直線（曲線）を描く。
４）式を表示する。
５）二乗誤差を計算して表示する。
プログラム
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」になっている。
図１ 得られた結果の一例
「Visual C++ の勉強部屋」（目次）へ