Apache Commons Math3 (3.4)を使ったLevenbergMarquardt (LM)法によるCurve fittingの方法 (1変数関数)

## はじめに

ネット上にそのものズバリの例が載ってないし、[Appach Commons Mathのページ](http://commons.apache.org/proper/commons-math/userguide/fitting.html)を見ても、いまいちよくわからないのでなんとか自分で作ってみる。

FijiのExample source codeに[Optimization_Example.java](https://code.google.com/p/fiji-bi/source/browse/src-plugins/Script_Editor/templates/Java/Optimization_Example.java?r=d0213fe689866e101ec55b092b4ef9d5276a5307)があるが、古い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`という簡単な関数を定義している。

```java
    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がない)である。このインターフェースは二つのメソッドを定義している。

```java
    double[] gradient(double x, double... parameters)
    double value(double x, double... parameters)
```    
gradientは「勾配」である。このメソッドでは、各パラメータで偏微分した値を、double型配列に格納する形で返す。プログラム中では次のように例として記述されている。

```java
    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 でもびっくりマークでてるし)。

```java
    public double[] gradient(double x, double... params) {
```

valueはその関数そのものを指す (引数は新しい形になおした)。

```java    
    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を初期化する必要があるらしい。

```java
    59: // Now we initialize the optimizer and curve fitter.
```

ここでは、

```java
    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**になっていて、次のように記述されている。

```java
    Deprecated. 
    As of 3.3. Please use AbstractCurveFitter and WeightedObservedPoints instead.
```
    
これをみると、AbstractなクラスであるAbstractCurveFitterを継承するクラスを作る必要があるらしい。もう一つの、WeightedObservedPointsは、データを登録する際に使われる。CurveFitterでは、データは データは次のように追加していた。

```java
    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](http://commons.apache.org/proper/commons-math/userguide/fitting.html)で説明されている。

```java
    // 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を修正すると、

```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](http://commons.apache.org/proper/commons-math/apidocs/org/apache/commons/math3/fitting/SimpleCurveFitter.html)というのを見つけた。このクラスの説明をみると、

```
    Fits points to a user-defined function.
```
    
とあるので、これを使えば良いみたい。SimpleCurveFitterにはクラスメソッドとして、`create(ParametricUnivariateFunction f, double[] start)`がある。引数の第１項は、ParametricUnivariateFunctionであり、定義した関数をここに入れればパラメータ化された一変数関数を作ってくれるみたい。次の引数のdouble型配列は初期値だ。

以上をまとめると、Optimization_Example.javaでは以下のようにすればよいはずだ。

```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](https://github.com/apache/commons-math/blob/master/src/test/java/org/apache/commons/math3/fitting/SimpleCurveFitterTest.java)がありました。でも、このプログラムの中身をみても、多分私には理解できていないだろう。この例では、PolynominalFunctionを使ってパラメトリックな関数を作って、SimpleCurveFitterに放り込んでいるだけなので、自分の好きな関数を使いたい時はどうしたらいいのかきっと悩んでたと思う。

次は多変数の場合のフィッティング方法を調べて、２次元ガウス関数のフィッティングができるようにする。

### 疑問点
今のところ以下の点がわからない。

- フィッティングした時のIterationの回数を知る方法はないか？（どれくらいで収束したか知りたい）。
- LM以外のoptimizerを指定する方法はあるか？
- ParametricUnivariateFunctionからRealの表記が消えた理由は何か？

・・・まあおいおいわかるだろう。

## 追記
Fijiの[Optimization_Example.java](https://github.com/fiji/Script_Editor/blob/master/src/main/resources/templates/Java/Optimization_Example.java)の、Apache Commons Math3に対応したコードがアップされていた(
気づかなかった・・・）。2014.Octに公開されているようだが、`org.apache.commons.math3.optimization.fitting.CurveFitter`をまだ使っているようだ。でも、そのうち新しいライブラリに対応したバージョンがアップされるだろう。

## 参考にしたサイト
- [Optimization_Example.java](https://code.google.com/p/fiji-bi/source/browse/src-plugins/Script_Editor/templates/Java/Optimization_Example.java?r=d0213fe689866e101ec55b092b4ef9d5276a5307)
- [17 Curve Fitting](http://commons.apache.org/proper/commons-math/userguide/fitting.html)
- [Apache Commons Math 3.4 API](http://commons.apache.org/proper/commons-math/apidocs/)
- [SimpleCurverFitterTest.java](https://github.com/apache/commons-math/blob/master/src/test/java/org/apache/commons/math3/fitting/SimpleCurveFitterTest.java)