最適化セミナーのご案内

2.16 最小二乗問題

 最小二乗問題は,科学や工学の分野でよく利用される基本的な問題です.この問題の典型的な例としては,曲線の当てはめがあります.実験により得られたデータをモデル曲線に当てはめ定数を推定していくことを考えます.すると,一般に各観測点においてデータとモデルの間には誤差が発生してしまいます.この時,「各観測点における誤差の二乗和を最小にする」という基準を採用しなるべく良い曲線に当てはめようとすると最小二乗問題となります.

例題

 Aさんは,実験開始からt秒後のyという値を計測するという実験を行いました.その結果をまとめると,下の表のようになりました.Aさんはこの結果をプロットしたところ,$f(t)=a{{t}^{2}}+bt+c$という二次関数モデルに当てはめられることに気付きました.そこで,各計測時刻でのデータ値とモデルから得られる値との誤差の二乗和が最小になるようにa, b, cを推定してください.

t y
0.5 1.22470
1.0 0.87334
1.5 0.99577
2.0 1.34215
2.5 2.12172

 

t y
3.0 3.22933
3.5 4.44744
4.0 6.13509
4.5 8.14697
5.0 10.50759

 この例題において,変数となるのは推定するa, b, cの3つです.また,ある計測時刻$t$でのデータ値$y(t)$とモデルから得られる値$f(t)$との誤差は$at^2 + bt + c - y(t)$となります.よって,各計測時刻について誤差の二乗を求め,その総和をとることで最小化する目的関数が得られます.

 以上のことから,この例題は次のように定式化できます.

変数
a 2次の項の係数
b 1次の項の係数
c 定数項
 
目的関数(最小化)
$(0.5^2a + 0.5b + c - 1.22470)^2 + (1.0^2a + 1.0b + c - 0.87334)^2 + (1.5^2a + 1.5b + c - 0.99577)^2 + (2.0^2a + 2.0b + c - 1.34215)^2 + (2.5^2a + 2.5b + c - 2.12172)^2 + (3.0^2a + 3.0b + c - 3.22933)^2 + (3.5^2a + 3.5b + c - 4.44744)^2 + (4.0^2a + 4.0b + c - 6.13509)^2 + (4.5^2a + 4.5b + c - 8.14697)^2 + (5.0^2a + 5.0b + c - 10.50759)^2$ 各計測時刻における誤差の二乗和

 ここで,この結果をSIMPLEで何も工夫せずに記述すると次のようになります.

//変数の宣言
Variable a(name="2次の項の係数");
Variable b(name="1次の項の係数");
Variable c(name="定数項");
//誤差の二乗和の最小化
Objective err(type=minimize);
err=(0.5*0.5*a+0.5*b+c-1.2247)*(0.5*0.5*a+0.5*b+c-1.2247)
+(1.0*1.0*a+1.0*b+c-0.87334)*(1.0*1.0*a+1.0*b+c-0.87334)
+(1.5*1.5*a+1.5*b+c-0.99577)*(1.5*1.5*a+1.5*b+c-0.99577)
+(2.0*2.0*a+2.0*b+c-1.34215)*(2.0*2.0*a+2.0*b+c-1.34215)
+(2.5*2.5*a+2.5*b+c-2.12172)*(2.5*2.5*a+2.5*b+c-2.12172)
+(3.0*3.0*a+3.0*b+c-3.22933)*(3.0*3.0*a+3.0*b+c-3.22933)
+(3.5*3.5*a+3.5*b+c-4.44744)*(3.5*3.5*a+3.5*b+c-4.44744)
+(4.0*4.0*a+4.0*b+c-6.13509)*(4.0*4.0*a+4.0*b+c-6.13509)
+(4.5*4.5*a+4.5*b+c-8.14697)*(4.5*4.5*a+4.5*b+c-8.14697)
+(5.0*5.0*a+5.0*b+c-10.50759)*(5.0*5.0*a+5.0*b+c-10.50759);
//求解し結果を表示する
solve();
a.val.print();
b.val.print();
c.val.print();

 これでは,何回も実験を行うような場合モデルファイルの修正に大変な手間がかかってしまい効率的ではありません.そこで,SIMPLEの特長を生かし様々な測定結果に対応できるように定式化から見直していきます.

 まず,ある計測時刻tでのデータ値$y(t)$とモデルから得られる値$f(t)$との誤差を$e(t)$と表すことにします.すると,先ほどの議論から$e(t) = at^2 + bt + c - y(t)$となります.これより,目的関数は$e(t)^2$をtについて和をとったものと言えます.SIMPLEにおいては,Expressionを用いることで数式を宣言できます.よって,$e(t)$をExpressionで記述することでモデルファイルをより簡潔なものにできます.また,二乗和を求める部分についても観測時間の集合を導入することでsum()を用い簡単に記述できます.最後に,各計測時刻tにおける$y(t)$の値は直接モデルファイルに記述するのではなくcsvファイルから与えることにします.

 以上のことから,次のように定式化しなおすことができます.

集合
$ObserveTime = \{ 0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,\allowbreak 4.5,5.0 \}$ 観測時間の集合
 
変数
a 2次の項の係数
b 1次の項の係数
c 定数項
 
定数
$y(t),t \in ObserveTime$ 時刻$t$における観測値
 
目的関数(最小化)
$\displaystyle \sum_{t \in ObserveTime} e(t)^2 = \sum_{t \in ObserveTime} (at^2 + bt + c - y(t))^2$ 各観測時刻における誤差の二乗和

 この結果をSIMPLEで記述すると次のようになり,先ほどのものに比べ汎用性が高くさらに見やすいものになっています.

//観測時間の集合の宣言
Set ObserveTime;
Element t(set=ObserveTime);
//観測値をパラメータとして宣言
Parameter y(index=t);
//変数の宣言
Variable a(name="2次の項の係数");
Variable b(name="1次の項の係数");
Variable c(name="定数項");
//誤差をExpressionを用いて表現する
Expression e(index=t);
e[t]=a*t*t+b*t+c-y[t];
//誤差の二乗和の最小化
Objective err(type=minimize);
err=sum(e[t]*e[t],t);
//求解して結果を出力する
solve();
a.val.print();
b.val.print();
c.val.print();

 なお,実行させる際には次のようなcsvファイルを与えます.

t,y
0.5,1.22470
1.0,0.87334
1.5,0.99577
2.0,1.34215
2.5,2.12172
3.0,3.22933
3.5,4.44744
4.0,6.13509
4.5,8.14697
5.0,10.50759

 このモデルを実行すると,最後に

2次の項の係数=0.644402
1次の項の係数=-1.47656
定数項=1.76057

という表示がされます.このことから,$f(t) = 0.644402t^2 - 1.47656t + 1.76057$という二次関数モデルが推定できたことになります.


 

 

上に戻る