はじめに
ネット上にそのものズバリの例が載ってないし、Appach Commons Mathのページを見ても、いまいちよくわからないのでなんとか自分で作ってみる。
FijiのExample source codeにOptimization_Example.javaがあるが、古いApache Commons mathを使っており、将来的に使えなくなるメソッドを使っている。したがって、math3に対応した形に書き直す必要がある。
Optimization_Example.javaでは、public void fit()とpublic void optimize()という二つの関数がある。fit()では一変数関数 (univariate function)をLM法でfittingする、ということを行っている。optimize()の方は２変数での関数のフィッティング例を示しているようだ。このうちfit()の方を眺めていく。
1変数関数によるCurveFitting
関数の定義
プログラム中のコメントによると、まず、パラメータを持つ一変数関数を定義する必要が有り、さらにその関数は微分可能である必要があるとのことである。ちなみに関数の例として、a*log(x)+bという簡単な関数を定義している。
36: /*
37: * First, we define a parameterized univariate function.
38: * To use the Levenberg-Marquardt strategy, it has to be differentiable.
39: * For demonstration purposes, we define it as a simple logarithmic function:
40: * f(x) = a * log(x) + b
41: */
その関数はインターフェースとしてorg.apache.commons.math.analysis.ParametricUnivariateRealFunctionを実装したクラスで行っているが、これは古いインターフェースである。(多分)、新しいインターフェースはorg.apache.commons.math3.analysis.ParametricUnivariateFunction(Realがない)である。このインターフェースは二つのメソッドを定義している。
double[] gradient(double x, double... parameters)
double value(double x, double... parameters)
gradientは「勾配」である。このメソッドでは、各パラメータで偏微分した値を、double型配列に格納する形で返す。プログラム中では次のように例として記述されている。
43: public double[] gradient(double x, double[] params) {
44: double a = params[0];
45: double b = params[1];
46: return new double[] {
47: Math.log(x),
48: 1
49: };
50: }
47行目が関数a*log(x)+bをaで偏微分した値、48行目がbで偏微分した値になっている。
ちなみに新しいライブラリでのパラメータは可変長引数なので、以下のように書き直すべきだろう (Eclipse でもびっくりマークでてるし)。
public double[] gradient(double x, double... params) {
valueはその関数そのものを指す (引数は新しい形になおした)。
52: public double value(double x, double... params) {
53: double a = params[0];
54: double b = params[1];
55: return a * Math.log(x) + b;
56: }
OptimizerとCurve Fitterの初期化??
つぎに、optimizerとcurve fitterを初期化する必要があるらしい。
59: // Now we initialize the optimizer and curve fitter.
ここでは、
60: LevenbergMarquardtOptimizer optimizer = new LevenbergMarquardtOptimizer();
61: CurveFitter fitter = new CurveFitter(optimizer);
としている。optimizerでは、どの方法でデータに対してモデルを最適化するのかを選び、curve fitterで具体的なfitting操作をするようだ。ここで、LevenbergMarquardtの文字が出てくる。古いライブラリでは、org.apache.commons.math.optimization.general.LevenbergMarquardtOptimizerを使っているが、現在のライブラリではorg.apache.commons.math3.fitting.leastsquares.LevenbergMarquardtOptimizerを使うようだ(
同じ名前のメソッドがライブラリ中に３つもあり、そのうち二つがDeprecatedとしてVersion4では消えるようだ)。
さて、問題は次のCurveFitterである。CurveFitterはライブラリ中に二つあるが(~.fittingパッケージににある)、どちらもDeprecatedになっていて、次のように記述されている。
Deprecated.
As of 3.3. Please use AbstractCurveFitter and WeightedObservedPoints instead.
これをみると、AbstractなクラスであるAbstractCurveFitterを継承するクラスを作る必要があるらしい。もう一つの、WeightedObservedPointsは、データを登録する際に使われる。CurveFitterでは、データは データは次のように追加していた。
63: // and feed it some "observed" data points of the form (x, y)
64: double[] x = { 1.1, 20.2, 100.3 };
65: double[] y = { 5.9, 4.8, 3.7 };
66: for (int i = 0; i < x.length; i++)
67: fitter.addObservedPoint(x[i], y[i]);
WeightedObservedPointsによるデータ登録
つまり、CurveFitterのメソッドであるaddObservedPoints(double x, double y)を使うことで、データを登録していた。新しいライブラリでは、WeightedObservedPointsを次のように使うことでデータを登録すべきと、Apache Commons Mathの17 Curve Fittingで説明されている。
// Collect data.
final WeightedObservedPoints obs = new WeightedObservedPoints();
obs.add(-1.00, 2.021170021833143);
obs.add(-0.99, 2.221135431136975);
obs.add(-0.98, 2.09985277659314);
obs.add(-0.97, 2.0211192647627025);
// ... Lots of lines omitted ...
obs.addt(0.99, -2.4345814727089854);
// Instantiate a third-degree polynomial fitter.
final PolynomialCurveFitter fitter = PolynomialCurveFitter.create(3);
// Retrieve fitted parameters (coefficients of the polynomial function).
final double[] coeff = fitter.fit(obs.toList());
最初にfinalとしてWeightedObservedPointsのオブジェクトobsを作る。そして、addメソッドを使ってデータを追加していく。このデータ登録例にならってOptimization_Example.javaを修正すると、
final WeightedObservedPoints obs = new WeightedObservedPoints();
double[] x = { 1.1, 20.2, 100.3 };
double[] y = { 5.9, 4.8, 3.7 };
for (int i = 0; i < x.length; i++)
obs.add(x[i],y[i]);
となる。
次に、CurveFitterとしてPolynominalCurveFitterを使い、３次の多項式を
createメソッドで作成している。最後に、fitメソッドでフィッティングを行う。引数としてWeightedObservedPointsのtoListメソッドの返り値を渡すことで、double型の配列として係数を得ることができるようだ。CurveFitterでは、多項式関数用のPolynominalFitter, 調和関数用のHarmonicFitter、そしてガウス関数用のGaussianFitterが標準で用意されている。確かにこの辺りがよく使われるのだろう (Exponentialくらいはあっても良い気がするが)。それ以外の関数はAbstractCurveFitterを使って定義する必要があるようだ。
SimpleCurveFitterによるユーザー定義関数によるフィッティング
・・・と思ったら、org.apache.commons.math3.fitting.SimpleCurveFitterというのを見つけた。このクラスの説明をみると、
Fits points to a user-defined function.
とあるので、これを使えば良いみたい。SimpleCurveFitterにはクラスメソッドとして、create(ParametricUnivariateFunction f, double[] start)がある。引数の第１項は、ParametricUnivariateFunctionであり、定義した関数をここに入れればパラメータ化された一変数関数を作ってくれるみたい。次の引数のdouble型配列は初期値だ。
以上をまとめると、Optimization_Example.javaでは以下のようにすればよいはずだ。
SimpleCurveFitter scf = SimpleCurveFitter.create(function, new double[] {1,0});
double[] result = scf.fit(obs.toList());
あれ、こんだけ？optimizerは指定しなくてよいのか？AbstractCurveFitterのgetOptimizer()の項目を見ると、次のように記述されている。
The default implementation uses a Levenberg-Marquardt optimizer.
つまり、デフォルトでLM法が選択されているので、LM法を使う限りにおいて、特にoptimizerを指定する必要はないらしい(逆に、optimizerを別のものにする場合はどうするんだろう？)。
まとめ
以上をまとめると、
-
org.apache.commons.math3.analysis.ParametricUnivariateFunctionを使って、関数を定義する -
org.apache.commons.math3.fitting.WeightedObservedPointで生データを登録する -
org.apache.commons.math3.fitting.SimpleCurveFitterのクラスメソッドcreateを使って、作った関数と初期値を用いたSimpleCurveFitterオブジェクトを作る。 -
SimpleCurveFitterのfitメソッドに、WeightedObservedPointで登録した生データをtoListメソッドをつかって渡してやるとフィッティングを行う。返り値にフィッティングしたパラメータが入る。
・・・めちゃ簡単や。abstractCurveFitingの記述に惑わされてしまった感がある。
ちなみに、「もしかして」と"apache commons math simplecurvefitter"でググると、SimpleCurverFitterTest.javaがありました。でも、このプログラムの中身をみても、多分私には理解できていないだろう。この例では、PolynominalFunctionを使ってパラメトリックな関数を作って、SimpleCurveFitterに放り込んでいるだけなので、自分の好きな関数を使いたい時はどうしたらいいのかきっと悩んでたと思う。
次は多変数の場合のフィッティング方法を調べて、２次元ガウス関数のフィッティングができるようにする。
疑問点
今のところ以下の点がわからない。
- フィッティングした時のIterationの回数を知る方法はないか？（どれくらいで収束したか知りたい）。
- LM以外のoptimizerを指定する方法はあるか？
- ParametricUnivariateFunctionからRealの表記が消えた理由は何か？
・・・まあおいおいわかるだろう。
追記
FijiのOptimization_Example.javaの、Apache Commons Math3に対応したコードがアップされていた(
気づかなかった・・・）。2014.Octに公開されているようだが、org.apache.commons.math3.optimization.fitting.CurveFitterをまだ使っているようだ。でも、そのうち新しいライブラリに対応したバージョンがアップされるだろう。