Apache Commons Mathを使って最小二乗法を試した

Apache Commonsは有名だけどなぜかMathのサンプルは見つからないんですね。。。なぜだろう？

とりあえず使ってみました。何かの役には立つでしょう。たぶん。。。

絵的なものが必要なのでProcessingで動かしました。

コードは以下におきました。よければ参考にしてください。

実行結果は以下のような感じです（と言われてもわかりにくいですね。。。最初は真っ暗なんですが、左クリックで点を追加していき右クリックを押すと当てはめ計算をします。計算結果でy=ax+bの係数a,bが得られますので、その値にしたがってlineで直線を描画しています。

f:id:ke_takahashi:20160314175932p:plain

以上です。

あ、少し補足しますね。

最小二乗法はもちろん知っていましたがなかなか使う機会がないので勉強もする機会もありませんでした。仕事でも学生のときも。

このコードはApache Commons Mathの導入ガイドのページのコードを改変しました。

ほとんど変えていないといってもいいと思います。このページにコードとライブラリの使用方法について書かれています。じっくり読むとどこを変えればよいか見えてきます（これだけではわからないので最小二乗法の勉強を少ししなければなりませんが）。

あとはこれを読んでね。。。ではあれなので（未来の僕も忘れているにちがいないので）変更したポイントを書き残したいと思います。

まず、observedPointsがそのままの訳で観測値ですね。フィッティングしたい値です。制限はないようなので何個でもいいんじゃないかと。

次にちょっと飛ばしてprescribedDistancesですね。これは観測値とフィッティングした結果求められるモデル式との差異です。この例では直線に当てはめたいので０にします。観測値それぞれに必要なので、観測値の個数分、領域を確保してそこに値をセットする必要があります。

で、その下の.start(new double[] { 0.0, 0.0 })ですね。これは探索するときの開始値です。今回はaとbの２つの係数の値を求めたのですが、それぞれ０の値でスタートさせています。

ここまではコードと解説ページをみれば何となくわかったのですが、MultivariateJacobianFunctionクラスのvalueメソッドが？でした。このメソッドは最適化計算中に呼び出されるものです。valueメソッドの引数であるpointにベクトル形式で計算途中の値（今回の場合のaとbが該当します）が格納されています。

（GAでいうところの評価関数になると思うのですが）今回のpointは最適値と比べてどれくらい良いかを求める必要があります。具体的にはmodelIとjacobianです。modelIには最小化問題の評価式を、jacobianには求めたい係数で偏微分した値を入れます。言葉の使い方が怪しいですが。。。数式で書きますね。

${\displaystyle modelI = \frac{1}{2} \sum_{k=1}^{N} (y_k-(ax_k+b))^{2} }$
${\displaystyle \frac{\partial modelI}{\partial a}=a \sum_{k=1}^{N}x_k^{2}+b \sum_{k=1}^{N}x_k- \sum_{k=1}^{N}x_ky_k }$
${\displaystyle \frac{\partial modelI}{\partial b}=a \sum_{k=1}^{N}x_k+b \sum_{k=1}^{N}1- \sum_{k=1}^{N}y_k }$