理科系の勉強日記

Linux/Ubuntu/Mac/Emacs/Computer vision/Robotics

統計的手法

数学
最小二乗法

n個の観測値の組{\mathbf{x}_i, y_i }が与えられているとき，2つの変量間の関係を説明するモデルとして


{ y = f( \mathbf{x}, \mathbf{\theta} ) + \epsilon }
を考える．ここで，{\mathbf{\theta}}はモデルのパラメータである．

平均2乗誤差(Least Mean Squared Error)を最小とするようなパラメータを求める推定方法は，最小2乗法と呼ばれている．つまり，


{LMS = \frac{1}{n} \sum_{i=1}^n \epsilon^2 = \frac{1}{n} \sum_{i=1}^n |y - f( \mathbf{x}, \mathbf{\theta} )|^2 }
である．

M-estimator

当てはめるモデルと観測値との誤差が平均で0の正規分布に従う場合，最小2乗法で推定したモデルは最適となる．しかし，観測値に外れ値（例外値）が含まれている場合，推定結果に大きな影響を与えてしまう．

そこでロバストな推定方法を用いることで，外れ値を多く含む観測値に対しても比較的良い推定結果を得ることができる．代表的な手法としてM-estimatorとLeas Median of Squaresなどがある．

最小2乗法で用いられる最小2乗基準


{LMS = min \sum_i \epsilon_i^2 }
の代わりに，少しだけ変形した評価基準を用いる．この手法をM-estimatorという．この場合，すべてのデータに対して均等に重みをつけているため，推定結果が外れ値に大きな影響を受けてしまう．例外値に小さな重みを与えるような評価基準として

{M = min \sum_i \rho(x) }
を扱う．{\rho(x)}は，0を唯一の最小値とする偶関数である．

{\rho(x) = x^2 } とした場合が最小2乗法である．つまり，M-estimatorは最小2乗法の拡張と言える．

{ \Psi(x) = \frac{ {\mathrm d} \rho(x)} { {\mathrm d} x } }はinfluence functionと呼ばれ，観測値がモデルから離れた場合の重みを表す．図1にその一例を示す．{\rho}関数のとりかたによっては，観測値がモデルからある程度離れるとその影響はほぼ0となることがわかる．


図1: (左上) {\rho(x) = x^2}， (右上){\Psi(x) = 2x }，(左下) {\rho(x) = \frac{x^x}{ \sigma + x^2}}， (右下) {\Psi(x) = \frac{2x\sigma}{(\sigma + x^2)^2}}

M-estimatorによる推定のためのアルゴリズムはMを最小化する最適化問題として定式化することができ，重み付き最小2乗問題となる．ただし，このアルゴリズムは，{\rho}関数によっては必ずしも最適解に収束するとは限らないため，良い初期値から出発する必要がある．

シンボリック変数

matlab

matlabを用いて，解析的に積分をする上で必要となったのでメモ．

シンボリック変数を用いて数式を定義しておけば，matlabを使って解析的な処理が可能となる．（曖昧）
まず，xの関数yを定義するために，以下のコマンドを入力．

syms x y
y = sin(x)^2

例として，関数yをxで微積分してみる．

diff(y,x)
 
ans =
 
2*cos(x)*sin(x)
int(y,x)
 
ans =
 
x/2 - sin(2*x)/4
int(y,x,0,pi)
 
ans =
 
pi/2

上から順に微分，不定積分（原始関数），定積分となっていることを確認．

これらをグラフ上にプロットする場合は

ezplot(y)
ezplot(diff(y,x))
ezplot(int(x,y))

などとすれば良い．

積分結果が初等関数で表せない積分

matlab

定積分


{ \int_0^a x^2 e^{-x^2} dx}
を計算しようとして詰まった．解析的に定積分を行うのは久々である．
"expの微分は，指数部分の微分が係数になる"という記憶しかない．

ここで，すべての数式が初等関数で表現できるとは限らないということを思い出す．定積分は存在するが，初等関数で積分した後の関数を表せない積分も存在する．

初等関数（しょとうかんすう）とは，複素数を変数とする多項式関数・指数関数・対数関数主値の四則演算・合成によって表示できる関数である．（Wikipedia 初等関数）
初等関数の導関数は必ず初等関数になるが、初等関数の原始関数、及び初等関数を用いた微分方程式の解は必ずしも初等関数になるとは限らない。 (Wikipedia 初等関数)

積分結果が初等関数で表せない場合，

    1. 数値積分を行う
    2. 特殊関数を使って積分結果を表す

という方法がある．

先に上げた数式の不定積分を後者の方法で表現すると以下のようになる．


{ \int_0^a x^2 e^{-x^2} dx = \frac{ \sqrt{\pi} } {4} erf( x ) - \frac{x}{2} e ^{-x^2} + C }
ただしCは定数，erfは誤差関数と呼ばれる特殊関数であり，

{ erf( x ) = \frac{2}{\pi} \int_0^x e^{-t^2} dt }
で表される．

matlabでplotするときに次元を縮退させる話 

hoge        <300x300x256>

というデータの3列目のグラフの概形が見たい．
前2つの最初の成分に対する3列目の値を二次元プロットする．

plot( hoge(1, 1, : ) )

Error using plot
Data may not have more than 2 dimensions

次元数が多いのでダメらしい．納得がいかないが仕方ないのでググってみて解法を得た．

plot( sqeeze(hoge(1, 1, : ) ) )

sqeezeは，引数に指定したベクトルの次元数を下げるコマンドである．
これで無事に列ベクトルになってくれたのでplotすることができた．

matlabエラーメッセージの怪

matlab

エラーメッセージに騙されて格闘した30分間に記録．


正規分布をつくるために，横軸0から255を256段階で用意し，標準偏差を適当にきめ，平均をグレースケールの画像上のある画素の輝度値から決めた．
normpdfは正規分布（matlab的には正規確率密度関数）を返す関数である．

img_gray = imread('sample.png');

lin = linspace(0,255,256);
sigma = 10;
mu = img_gray(1,1);

pdfdata = normpdf(lin, mu, sigma);

すると以下の様なエラーが返ってきた．

Error using normpdf (line 36)
Non-scalar arguments must match in size.

normpdfを使う上でエラーが発生した．スカラじゃない引数はサイズが一致している必要があるという．
もう一度引数を確認すると

lin <1x256 double>
mu 16
sigma 10

となっており，問題があるように思えない．試しに以下のようにmuを宣言して実行してみる．

mu = 16;
pdfdata = normpdf(lin, mu, sigma);

するとうまくいった．直接値をmuに代入するとうまくいくようだ．
おそらくスカラとベクトルの問題ではなく，データ型の問題のような気がしてくる．
エラーメッセージより，normpdf.mの36行目付近を見にいく．

try
    y = exp(-0.5 * ((x - mu)./sigma).^2) ./ (sqrt(2*pi) .* sigma);
catch
    error(message('stats:normpdf:InputSizeMismatch'));
end

この条件分岐では，エラーが出たときはいつも「サイズのミスマッチ」と表示するようだ！
上の数式だけをコピーして実行すると

Undefined function 'exp' for input arguments of type 'uint8'.

というエラーメッセージが返ってきた．やはりデータ型の問題だったようだ．
今回グレースケールの画像は符号無し整数型で読み込んだが，expはdoubleを引数にとるんだね．

mu_double = cast(mu, 'double)'; 
pdfdata = normpdf(lin, mu_double, sigma);

これで上手くいった！
エラーメッセージを信用してはならないという教訓を得た．