2017-01-31

MATLABで非線形最小二乗フィッティングする手順メモ

多数のデータ群から関数フィッティングを行う場合に非線形最小二乗法を用います。

MATLABにはOptimization Toolboxなどがあり，それを用いると簡単にフィッティングできます。
以下のブログなんかに解説がありますね。
MATLABで非線形最小二乗問題を解く - Py

このようなToolboxなしでも，
地道にこれを解く手順がありますので勘所をメモしておきます。
なお，数式とアルゴリズムの基礎がわかっている前提で書きます。

ダンピングを入れたVerはこちら。
ossyaritoori.hatenablog.com

アウトライン
問題設定
問題の定式化
ヤコビ行列の作成
- 補足：ヤコビ行列のサイズ
- 注意：subsを用いた複数の配列の代入
反復計算による求解
- 補足：行列のサイズ
計算結果
まとめと今後
サンプルコード
参考資料
jacobian 関数の使用：2017/3/7 更新

アウトライン

matlabのバージョンは生協で買ったStudent version（2015a準拠)です。
本記事ではsymbolic math Toolboxを使って「偏微分」と「変数代入」を行います。
（頑張ればここは手作業でいけます。しないけど。）

ヤコビアンを求めて，Gauss-Newton法の枠組みでこれを計算します。

こんな風にフィッティング出来ます。
f:id:ossyaritoori:20170131220211p:plain

なお，Gauss-Newton法は案外簡単に発散するのでその場合は以下のLevenberg-Marquardt法を試してみてください。
ossyaritoori.hatenablog.com

問題設定

例題としてwikipedeiaの問題を考えます。
ガウス・ニュートン法 - Wikipedia

Fitting対象の関数は以下のようになります。

$y = \dfrac{ax}{x+b}$

フィッティングするのはaとbの2つの変数です。

フィッティングするデータは次の様に表せます。

% data for fitting
xhat = [0.038, 0.194, 0.425, 0.626, 1.253, 2.500, 3.740].';
yhat = [0.050, 0.127, 0.094, 0.2122, 0.2729, 0.2665, 0.3317].';

通常こういう値のベクトルは縦向きですので後ろで転置してます。
今回は7×1のベクトル，すなわち7個のデータセットですね。

問題の定式化

残差と最小二乗和を以下の様に定義します。iはデータ群の各成分ってこと。

残差： $r_i = y_i - \dfrac{ax_i}{x_i+b}$

残差の二乗和： $S = \sum_{i=1}^n (r_i)^2$

Symbolic math Toolboxを使ってこれは次の様にかけます。

%% making evaluate funciton
syms a b x y

% evaluate fuction
func = a*x/(x+b);
R = y - func;
S = R^2;

ヤコビ行列の作成

ヤコビ行列を求めます。偏微分はシンボリックの微分diffで簡単に作成できます。
（もっと簡単に作成する手法は一番下を参照）

実装上工夫する点は，

代入できる値は早めに代入して，計算コストを減らす，（シンボリック計算は重い）

% partial derivative
Ra = diff(R,a);
Rb = diff(R,b);

% make jacobian matrix
Ja = subs(Ra,{x,y},{xhat,yhat});
Jb = subs(Rb,{x,y},{xhat,yhat});

Jacob = [Ja,Jb]; % jacobian

Ja,Jbの計算では，ヤコビ行列に先にフィッティングデータの値を代入して後の計算を楽にします。

補足：ヤコビ行列のサイズ

フィッティング係数がm個，データがn個あるとします。

この際，ヤコビ行列はn×mになることに注意。
f:id:ossyaritoori:20170131214013p:plain:w200

注意：subsを用いた複数の配列の代入

subs関数を用いてシンボリックに配列を代入する時の書き方ですが，ややコツがいるかも？
公式を読んでいると以下のように代入しがちですがアウトです。

% 失敗例
subs(fa,[x,y],[xhat,yhat]);

セル配列として受け渡すのが正解です。

% 成功例
subs(fa,{x,y},{xhat,yhat});

セル配列誰かに教わりたい。

反復計算による求解

ヤコビ行列が求まれば，その時々のパラメータの推定値 $a_n,b_n$ をヤコビ行列Jと残差Rに代入して更新量がわかります。
パラメータをまとめてvとし，残差への代入値をrとすると，

更新式： $v_{new} = v + \Delta v$
に対して

$\Delta v = - (J^TJ)^{-1}J^T \ r$

が求める計算です。
これまたsubsでの代入に注意。

% 残差に事前にフィッティングするデータを代入。
r = subs(R,{x,y},{xhat,yhat});

%ヤコビ行列にデータvを代入
J = double(subs(Jacob,[a,b],v'));
% 残差にも代入
df = double(subs(r,[a,b],v'));

% 更新量Δを計算して更新
delta = - ( J.' * J )\J.' *  df;
v = v + delta;

symsだと計算した値がシンボリックのままになるので，
doubleで一度数値に直しています。

補足：行列のサイズ

フィッティング係数がm個，データがn個あるとします。
この例題では，m=2，n=7ですね。

図解すると次のとおりになります。
f:id:ossyaritoori:20170131215154p:plain

計算結果

今回は固定ステップ itaration=10 でやりました。
一般的には，「誤差が閾値を下回るまで」という条件付けをします。
正確さを気にしなければ誤差の微分見て鈍くなった所で止めてもいいかと。

フィッティングの結果を図示します。
f:id:ossyaritoori:20170131220211p:plain

各ステップ毎の誤差の二乗和Sは次のようになっています。
f:id:ossyaritoori:20170131220401p:plain

順調に減って居ますが途中から鈍くなってますね。今回最適化した結果は

a = 0.3618
b = 0.5563

です。wikiの値とほぼ同じですね。

まとめと今後

simbolicでGauss-newton法のフィッティングを行いました。
今後やるとしたら以下についてまとめる予定です。

Levenberg-Marquardt アルゴリズムの実装
有限差分ヤコビアンを用いる -> symbolic抜きに擬似(数値)微分

サンプルコード

終了条件が誤差の大きさではなく，反復回数になっているので気に食わない方は変更してね。

%% get data
clear all;
% data for fitting
xhat = [0.038, 0.194, 0.425, 0.626, 1.253, 2.500, 3.740].';
yhat = [0.050, 0.127, 0.094, 0.2122, 0.2729, 0.2665, 0.3317].';

%% making evaluate funciton
syms a b x y

% evaluate fuction
func = a*x/(x+b);
R = y - func;
S = R^2;

% partial derivative
Ra = diff(R,a);
Rb = diff(R,b);

% make jacobian matrix
Ja = subs(Ra,{x,y},{xhat,yhat});
% subs(fa,[x,y],[xhat,yhat]); ではダメ
Jb = subs(Rb,{x,y},{xhat,yhat});

Jacob = [Ja,Jb]; % jacobian


%% initial value
vinit = [0.9;0.2];
v = vinit;
iteration = 10;

% evaluate function
SumOfError = zeros(iteration,1);
SEF = subs(S,{x,y},{xhat,yhat});

%% iterative optimization
% error function
r = subs(R,{x,y},{xhat,yhat});

for i = 1:iteration
   J = double(subs(Jacob,[a,b],v'));
   df = double(subs(r,[a,b],v'));
   delta = - ( J.' * J )\J.' *  df;
   v = v + delta;

   SOE = double(subs(SEF,[a,b],v'));
   SumOfError(i) = SOE.' * SOE;
end

%% showing
fit = subs(func,x,xhat);
yfit = double( subs(fit,[a,b],v') );

figure(1);
plot(xhat,yhat,'+',xhat,yfit);
xlabel('x');
ylabel('y');
legend('data','fitted value','Location','Best');
grid on;

figure(2);
plot(SumOfError);
xlabel('iteration');
ylabel('Sum Of Squared Error');
legend('Squared Error value')
grid on;

参考資料

基本的には教科書を読もう。

Wikipediaの記事（わかりやすい）

ガウス・ニュートン法 - Wikipedia

pythonでWikipediaの練習問題を解いたブログ（有限差分ヤコビアンを使用。わかりやすい）

Python NumPy サンプルコード: Gauss-Newton 法で非線形最小二乗問題を解く | org-技術

jacobian 関数の使用：2017/3/7 更新

「ヤコビ行列はシンボリックの微分diffで簡単に作成できます。」
matlabにはjacobian関数が存在し，ヤコビ行列を一発で作成してくれます。

% partial derivative
Jacob = jacobian(R,[a b])；

Referenceをちゃんと読もうな。

2017-01-20

初歩過ぎてあまり説明を見ないベクトル・行列の微分の公式

行列微分勉強

行列で偏微分をしなきゃいけない時にとっさに忘れて非常に困ったので自分用にメモを残します。

ベクトルでの微分
行列での微分
- 公式
メモ

ベクトルでの微分

Aとxはベクトル（行列）とします。

１． $\dfrac{\partial A^T x}{\partial x} = A$

２． $\dfrac{\partial x^T A}{\partial x} = A$

手計算してチェックしてもいいですが，パっと思いだすのにはやっぱ書いておくべきでしょう。
もちろん，適当に成分を仮定して計算することも出来ます。

３． $\dfrac{\partial x^T A x}{\partial x} = (A+A^T) x$

４． $\frac{ \partial }{ \partial x } \left(Ax - b \right)^T C \left(Ax - b \right) = A^T \left(C+C^T\right)\left( A x - b \right)$

３の練習用の証明)
$\dfrac{\partial x^T A x}{\partial x} = \dfrac{\partial (x^T A) x}{\partial x}+\dfrac{\partial x^T (A x)}{\partial x}$

積の微分法。表記にはすこぶる自信がない。それぞれ１と２を適用して下のようになる。

$= (x^T A)^T + A x$

ここで， $(AB)^T = B^TA^T$ なので，

$= A^Tx + A x$

とまぁこんな感じか。

4も地道に展開すれば証明できます。

行列での微分

$\begin{equation} \dfrac{\partial}{\partial X} a^\top X b = ab^\top \end{equation}$

公式

以下のサイトより引用
「ベクトルで微分・行列で微分」公式まとめ - Qiita

${\begin{eqnarray*} \frac{ \partial }{ \partial A } tr \left( A \right) &=& I \\ \frac{ \partial }{ \partial A } tr \left( AB \right) &=& B^T \\ \frac{ \partial }{ \partial A } tr \left( BA \right) &=& B^T \\ \frac{ \partial }{ \partial A } tr \left( ABA^T \right) &=& A\left(B + B^T \right)\\ \frac{ \partial }{ \partial A } tr \left(f \left(U = g\left(A\right)\right)\right) &=& \frac{ \partial }{ \partial U } tr \left(f \left(U\right)\right) \frac{ \partial }{ \partial A } tr\left(g\left(A\right)\right) \\ \frac{ \partial }{ \partial A^T } f \left( A \right) &=& \left(\frac{ \partial }{ \partial A } f \left( A \right) \right)^T\\ \frac{ \partial }{ \partial A } |A| &=& |A|\left(A^{-1} \right)^T\\ \frac{ \partial }{ \partial A } \ln|A| &=& \left(A^{-1} \right)^T\\ \frac{ \partial }{ \partial X } ||AX - B||_F^2 &=& 2A^T\left(AX - B \right)\\ \frac{ \partial }{ \partial X } ||XA - B||_F^2 &=& 2\left(XA - B \right) A^T\\ \end{eqnarray*} }$

メモ

2ノルムの計算はそのまま行列の掛け算で分解するよりも内積の形に変更して分配・交換を駆使した方がきれいになる。

2017-01-15

MATLABでNaNやInfを何とかしたい

matlab

ダジャレみたいな軽い記事です。

NaNを0に置換
Infを0に置換
種明かし
- よくわからんので実験

NaNを0に置換

Aという行列の中にNanがあってそれを0にする場合，以下のように書きます。

A(isnan(A))=0

なお，matlabではfor文はべらぼうに時間がかかるので
なるべく行列のまま処理したほうが計算時間の節約になります。

Infを0に置換

Aという行列の中にInfがあってそれを0にする場合，以下のように書きます。
なお，-Infも同様です

A(~isfinite(A))=0

種明かし

isnan(A)

というコマンドは行列Aのうち，Nanがある所を1，無い所を0の「論理値」で出力します。
はい，ここ重要です。
行列に論理値を食わせることで該当する箇所の値を参照できるらしいです。
なるほど～。

よくわからんので実験

以下のコードで実験してみました。

3*3の魔法陣を作ります。

Magi =[ 8  1  6; 3  5  7; 4  9  2]

代入するのは代数値でなければなりません。
例えば代数値の対角行列を入れてみる。

 Magi(diag([true true true]))

結果，ans = 8,5,2 と見事に対角の値だけにアクセスしています。

ならば，一つだけ食わせるとどうか？
なんと先頭の値にアクセスします。

 Magi(true)

結果：ans = 8

そして数値を足してずらせるようです。ポインタみたいだな。

 Magi(true+1)

結果：ans = 3

行 → 列の順に格納しているのがわかります。

 Magi(true+3)

結果：ans = 1

ほらね。

また，こんなことも出来ます。

Magi(true,true+1)

結果：ans = 1

~~う～ん，true=1として処理してるだけに見えてきた。~~

Infの方はなんか出力が逆ですが基本やってることは同じです。
ほ～ん。

2017-01-13

Tex ：定理や式分岐などを書くセット（dcasesを使おう）

Tex

よく忘れるのでこの際テンプレとしてメモ。
特にcasesで分岐のフォントが小さい時にdcasesというパッケージがある（＆styファイルが異なる）のは初耳でした。

定理を書くtheoremのテンプレ。screenの四角付き。

\begin{screen}
\begin{theorem}
	関数 $ f(x,y)$ において点 $ (a,b)$ が  $ f_x(a,b)=0$,  $ f_y(a,b)=0$ をみたすとき, $ f(a,b)$ が極値となるための判定条件は次の通りである.\\ 
	但し，$\displaystyle D(a,b)=f_{xy}(a,b)^2-f_{xx}(a,b)f_{yy}(a,b)$とする。
	\begin{enumerate}
		\renewcommand{\labelenumi}{(\roman{enumi}).}
	\item $ D(a,b)<0$,  $ f_{xx}(a,b)>0$ のとき, $ f(x,y)$ は点 $ (a,b)$ で極小値をとる.
	\item $ D(a,b)<0$,  $ f_{xx}(a,b)<0$ のとき, $ f(x,y)$ は点 $ (a,b)$ で極大値をとる.
	\item $ D(a,b)>0$ のとき, $ f(x,y)$ は点 $ (a,b)$ で極値をとらない.
	\item $ D(a,b)=0$ のとき,個別に判定する.
	\end{enumerate}
\end{theorem}
\end{screen}

こうなります。
f:id:ossyaritoori:20170113001621p:plain

式分岐を書く，caseのテンプレ。

\dfrac{\partial^2 \bm{J}}{\partial \theta_{1x}\partial \theta_{1y}} = 
\begin{cases}
-\omega_{yx} \leq 0 & (y < x, x \neq 2)\\
-\omega_{xy} \leq 0 & (y > x, x \neq n)
\end{cases}

こうなります。
f:id:ossyaritoori:20170113001702p:plain

場合分けでは見た目的にはdcasesを用いた方が良い

mathtool.styをインクルード。

\usepackage{mathtools}

dcasesと書くだけです。

\dfrac{\partial^2 \bm{J}}{\partial \theta_{1x}\partial \theta_{1y}} = 
\begin{dcases}
-\omega_{yx} \leq 0 & (y < x, x \neq 2)\\
-\omega_{xy} \leq 0 & (y > x, x \neq n)
\end{dcases}

この書き方だとわかりませんがｗ
Σなどを用いているなら違いがわかると思います。
f:id:ossyaritoori:20170113011020p:plain

2017-01-12

Tex 改行をまたいで括弧をつける

Tex

Tex小ネタです。

 \left\{ xxxxxxxxx \right\}

と書いて大きい括弧を書けるのは周知のとおりですが，これに対して次のように改行を挟むとエラーが出ます。

 \left\{ xxxxxxxxx \\
 yyyyyyyyy  \right\}

解決法は，\left.や\right.を用いて以下のように書きます。地味にめんどくせぇ…

 \left\{ xxxxxxxxx \right. \\
 \left. yyyyyyyyy  \right\}

2017-01-02

iTunesを用いたバックアップからの復元・データ移行 -セッションが失敗したため、iTunesはiPad”xxx”を復元できませんでした-の解決

バックアップ iTunes 復元自動クリック

実家に帰ったら親にiPadのデータ移行をしてほしいと頼まれました。
なんだそれだけかと思ったらこれがなかなかうまく行きません。

iTunesを用いたiPhone，iPadのバックアップ，復元の注意点
それでもうまくいかない！
『連打くん』を使った復元成功法
- コツ
- おまけ
まとめ

iTunesを用いたiPhone，iPadのバックアップ，復元の注意点

詳細な手順は他のサイトを見てもらうとして，この際に気をつけるべきは

- 使用しているiOS，iTunesは正しくアップグレードされているか
OSやiTunesのバージョンが異なるとエラーの元なので一度環境を整えてからやってみてください。

ちなみにソフトを最新にした後のウチの環境は，
PC：windows7
メモリ：２G
CPU：Corei3（4年前くらいの）
です。メモリ2Gって…orz

それでもうまくいかない！

うちはそれでもうまく行きませんでした。

復元をボタンを押してもしばらくすると
『iTunesはiPad”xxx”を復元できませんでした』のエラーが出ます。

そこで以下のサイトの情報を見つけました。
セッションが失敗したため、iTunesはiPhone”xxx”を復元できませんでした。iPhoneの復元ができない。 | Astro.jp

解決法は，，，

”バックアップから復元中”のwindowをクリック，ドラッグし続ける！

な，何を言っているのかわからないと思うが(以下略
f:id:ossyaritoori:20170102010912p:plain:w200

試しにwindowをポチポチクリックとドラッグしていたらこれが実際効いたんです！
対象のプロセスを常に最優先にしていないとリソースが足りなくなる的な何かなんでしょうか。

このままでは時間と労力をとても消費します…
うちのは1時間くらいかかったので人力ではやめたほうが良さそうです。

自動でクリックし続けるソフト（簡単）

クリック動作を自動でやってくれるものとして，
連打くんというツールがあるみたいです。
Vector：連打くん (WindowsMe/98/95用ソフト / ユーティリティ)

上記のURLでダウンロード可能です。

※フリーソフトなので自己責任で管理をお願いします。

簡単操作かつ強力だったので以降これでやっていきます。

f:id:ossyaritoori:20170102020324p:plain

UWSC を使う（経験者向け）

同種のソフトにUWSCがあります。こちらは任意にプログラム可能であるため，もっと凝った操作が可能です。
自動ツール作成経験者はこれを使ってもいいでしょう。（私も持ってますがやらない）

また，AutoClickPositionsというソフトを用いてみましたがクリック頻度が低くおすすめできません。

Pythonで書く（暇人向け）

Pythonをインスコしているなら自分で自動クリックプログラムを書くことができます。
以下に基本ルーティンを書きます。
マウスの座標をうまく取得してループで回せばきっとなんとかしてくれることでしょう。

正直pythonがインスコ出来てて，まともに動くPCなら上記問題は起きないかと思いますが。

# -*- coding: utf-8 -*-
from ctypes import *
user32 = windll.user32
import time
 
if __name__ == "__main__":
    #ここから処理
    x = 0 #座標指定
    y = 0
    user32.SetCursorPos(x,y) #カーソル座標移動
    user32.mouse_event(0x2,0,0,0,0) #左クリックを押す
    time.sleep(0.2) #0.2秒待機
    user32.mouse_event(0x4,0,0,0,0) #左クリックを放す

『連打くん』を使った復元成功法

予め言いますがこの作業は通しで10時間以上かかります。急ぎでなくかつPCを触らない就寝時などに行ってください。
連打くんを用いることで時間を大幅に短縮！
すごいぞ！連打くん！

プログラムを解凍してexeファイルを実効すると，
f:id:ossyaritoori:20170102020324p:plain
のような画面が出ます。

クリック間隔を1秒前後にして，「F9」を押すと自動でクリックし始めます。
従って入れる数字は1000より少し小さいくらいがおすすめです。
そして「F10」をおして終了します。

コツ

最初の「かかる時間を計算しています」のところが1番シビアな印象です。
ここは人力でwindowをクリックした方が良いかもしれません。

その後残り時間が表示されるようになったら連打くんにまかせる，といった手法を取りました。

おまけ

クリック間隔とかかる時間の関係。

クリック間隔：3秒 → かかる時間25時間（推定）
クリック間隔：0.8秒 → かかる時間1時間（推定）

なんなんだ。

まとめ

バックアップの復元で困ったときはウインドウをクリックし続けると良い。
フリーソフト「連打くん」を使って連打させよう。
古いパソコンは買い換えよう。
iTunesきらい。

2016-12-29

opencv_2410core.pdb not found とか出たら

安心してください。基本的にはあなたのプログラムの何処かでセグフォってます。
ただそれだけです。

このpdbファイルのサーバーとか関係なしにこのエラーがでることがあります。
何が言いたいかって言うとVisualStudioは滅べってこと。