Nuorium Optimizer C++SIMPLE例題集 V26

2.16 最小二乗問題

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

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

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

　ここで，この結果をC++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();

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

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

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

　この結果をC++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

という表示がされます．このことから，という二次関数モデルが推定できたことになります．

---

 

 

上に戻る