目的
Levenberg-Marquardt法を使用して非線形最小二乗問題を解き、近似曲線を求める。
使用ライブラリ
フリーで使用できるLevenberg-Marquardt法のC#ライブラリは、
のふたつがみつかった。
csmpfitはC言語で記述されたライブラリを文法だけC#に変更したもののようで、オブジェクト指向的なコーディングにはなっていないようだったが、ライセンスがMITライセンスであることが明確になっている点がありがたい。そのため、こちらを採用した。
Levenberg-Marquardt.NETはドキュメントが整備されている点が優れている。しかし、フリーであると謳われているものの、ライセンスが明記されておらず、また、VisualStudio2013で記述されているため、使用するVisualStudioのバージョンの問題で採用できなかった。
csmpfitの使用方法
テスト関数の流用
csmpfitを使うにあたり、困ったのは、ドキュメントがどこにあるかわからないことだ。
しかし、幸運にも、テストが同梱されている。
MPFitLib.TestプロジェクトのTestMPFit.csがそれだ。
ここには、
- 一次関数へのフィッティングをおこなうTestLinFit
- 二次関数へのフィッティングをおこなうTestQuadFit
- ガウス関数へのフィッティングをおこなうTestGaussFit
が記述されている。
また、
- 二次関数へのフィッティング時にパラメータを修正するTestQuadFix
- ガウス関数へのフィッティング時にパラメータを修正するTestGaussFix
というテストも用意されている。
これらを参考にすれば、ドキュメントなしでもどうにかできそうだ。
テスト関数以外に流用するコード
それから、このテスト関数を動かすため、テストプロジェクトにある、
- CustomUserVariable
- フィッティング対象となる関数
が必要になる。
前者は、CustomeUserVariable.csに記述されている。
後者は、ForwardModels.csに記述されているいくつかの関数から選ぶか、または、これらを参考に自分で記述することになる。
一次関数、二次関数については、そのまま流用できるだろう。
少し変更すれば、三次関数や一変数多項式が簡単に実現できる。
ただ、ガウス関数については、
y = N * exp(-0.5 * ((x - μ) ^ 2) / (σ ^ 2)) - a
のように、パラメータを4つ使用しているので、使用目的によってはこのままでは使えないかもしれない。
フィッティングの実行
テストコードを流用して、だいたい以下のようになった。
var v = new CustomUserVariable()
{
X = ..., // double[]型、測定値のXの値を格納
Y = ..., // double[]型、測定値のYの値を格納
Ey = ... // double[]型、値はすべて0.5,要素数は測定値の数だけ
};
var p = new double[] {
パラメータ0の初期値,
パラメータ1の初期値,
パラメータ2の初期値,
...};
var pars = new mp_par[p.Length];
for (int i = 0; i < pars.Length; i++)
pars[i] = new mp_par();
var result = new mp_result(p.Length);
int state = MPFit.Solve(
xxxxFunc, // 使用する関数
v.X.Length,
pars.Length,
p,
pars,
null,
v,
ref result);
if (state >= 1 && state <= 4)
// succeeded
return p;
else
return null;
まとめ
csmpfitを使用することで、C#によってカーブフィッティングをおこなうことができた。
テスト環境では、数百のデータを1ミリ秒かかるかどうか程度でフィッティングが終了した。
Levenberg-Marquardt法は、必ずしも毎回収束するわけではないが、それでも最小二乗法を解く優れたアルゴリズムであり、実験結果をまとめたい方には使用場面が多いように思う。
自分でガウス・ニュートン法などをプログラミングするよりは、こうしたライブラリを活用して楽をしたい。