目次
- 目次
- はじめに
- Ceres Solverによる曲線フィッティング
- Ceres Solverによるロバスト曲線フィッティング
- Ceres Solver関連記事
- 参考資料
- MyEnigma Supporters
はじめに
今回の記事もこれらの記事の続きです。
今回は、Ceres Solverを使って
曲線フィッティングをしてみようと思います。
Ceres Solverによる曲線フィッティング
先日までの記事で扱っていたサンプルでは、
データはなく、関数の最小値などを計算していました。
しかし、元々最小二乗法や非線形最適化が開発された目的は、
データの曲線当てはめでした。
今回は、このデータの曲線当てはめのサンプルコードを作成してみたいと思います。
今回は下記の指数関数に、
m=0,σ=0.2のガウス分布を付加したデータを
準備します。
今回は下記のコードを使って、
#!/usr/bin/env python # -*- coding: utf-8 -*- # author:Atsushi Sakai import matplotlib.pyplot as plt import numpy as np x=np.arange(0,10,0.1) y=[np.exp(0.3*tx+0.1)+np.random.normal(0.0,0.2) for tx in x] plt.plot(x,y,'+r') plt.grid(True) plt.show() f=open('data.csv', 'w') for (x,y) in zip(x,y): f.write(str(x)+","+str(y)+"\n") f.close()
下記のグラフのような、
データを生成しました。
続いて、このデータファイルを読み込ませて、
Ceresを使って曲線当てはめをしてみます。
まず初めにコスト関数を設定します。
それぞれのデータに対してコストを計算し、
その総和をコストとすると、下記のようになります。
/** * @brief コスト関数 **/ struct ExponentialResidual{ ExponentialResidual(double x, double y) :x_(x), y_(y){} template <typename T> bool operator()(const T* const m, const T* const c, T* residual) const { residual[0]=T(y_)-exp(m[0]*T(x_)+c[0]); return true; } private: //Observations for a sample const double x_; const double y_; };
データを使って、コスト関数を計算する場合は、
データを格納するメンバ変数を作り、
コンストラクタでその値を設定するような
コスト関数構造体を作る必要があります。
続いて、先ほどのコスト関数構造体を使って、
最適化の設定をします。
データのx,yのデータがそれぞれ格納された配列があるとして、
そのデータをそれぞれループで取得し、
先程のコスト関数に入力したオブジェクトを
AutoDiffCostFunctionでコスト関数オブジェクトに追加していく形です。
//最適化問題を解く変数と初期値の設定 double m = 0.0;//exp(mx+c) double c = 0.0; //最適化問題を解く用のオブジェクトの生成 Problem problem; for(int i=0; i<x.size(); i++){ CostFunction* cost_function=new AutoDiffCostFunction<ExponentialResidual, 1, 1, 1>( new ExponentialResidual(x[i],y[i])); problem.AddResidualBlock(cost_function, NULL, &m, &c); }
あとは以前の記事で紹介した最適化の方法と同じように、
最適化を実施すると下記のような結果が得られます。
上記の結果を見るとわかるように、
m=0.0, c=0.0から初めた最適化が、
m=0.300439, c=0.0966003と、
正解であるm=0.3,c=0.1にかなり近い値で計算出来ていることがわかります。
上記の最適化のコードはこちらで公開しています。
Ceres Solverによるロバスト曲線フィッティング
前述の方法は、
入力データがある程度正確であれば問題ありませんが、
もしデータの中にOutlierが含まれる場合(ノイズモデルに従わないデータ)、
下記の図のように、最適化された結果がOutlierに引きずられてしまいます。
そんな時に利用すると良いのが、損失関数です。
損失関数をうまく設定することで、outlierに強い最適化を実施することができます。
CeresSolverの場合、
AddResidualBlockの２つ目の引数を、
problem.AddResidualBlock(cost_function, NULL , &m, &c);
下記のように変更することにより、
problem.AddResidualBlock(cost_function, new CauchyLoss(0.5) , &m, &c);
損失関数を設定することができます。
CauchyLossはCeresソルバの損失関数の一つであり、
0.5は損失関数のスケールパラメータです。
ちなみにCauchyLossを使う場合は、
using ceres::CauchyLossを宣言するのを忘れないようにしましょう。
上記のように、損失関数を設定すると、
下記のようにOutlierが含まれるデータであっても、
真の値に近い、最適化を実施することができます。
上記の最適化のコードは下記で公開しています。
Ceres Solver関連記事
参考資料
MyEnigma Supporters
もしこの記事が参考になり、
ブログをサポートしたいと思われた方は、
こちらからよろしくお願いします。
