最適化計算２ 楕円の当てはめ 理論編
前回の下記投稿にひきつづき，今回は楕円の当てはめです．
daily-tech.hatenablog.com
１．楕円当てはめ問題
前回の投稿では，取得したデータをプロットし，そのプロットに対して直線を当てはめました．今回はデータに対して楕円を当てはめます．イメージとしては下記のような感じです．
前回の投稿では「一番それっぽい直線」を求めるために，最小二乗法を使いました．今回の楕円のケースでは，「一番それっぽい楕円」を求めるために下記の複数の方法を試します．
１．最小二乗法
２．最尤推定法
で，インプットとなる取得データとアウトプットとなる楕円の方程式は下記のようになります．
上式を見てみると，結局パラメータは増えているものの，直線当てはめの場合と似たような形の式になっていることがわかります．パラメータの正規化条件も加えて，下記のように整理します．
２．最小二乗法
最小二乗法を用いる場合考え方はシンプルで，直線当てはめの場合と同じ感じで式を立てます．つまり，取得データを楕円の方程式に代入し，この総和を最小にするようなパラメータを決めます．
上式は二次形式になっているので，直線当てはめの場合と同じように，実対称行列 MLS に対する最小固有値に対応する固有ベクトルが求めるパラメータベクトルになります．つまり，，，
ただーし，楕円の当てはめに対しては上記の最小二乗法では精度がそれほど出ません．下記エントリで最尤推定に関して扱いましたが，この時は直線の当てはめだったので "Ax + By + C = 0" というのが満たすべき式でした．ここで，"取得データ x, y が持つ誤差は互いに独立である"という過程のもと，尤度を最大化した結果が最小二乗法でした．
daily-tech.hatenablog.com
今回の楕円当てはめの場合では，x2，y2，xy，という項がありますが，この項に生じる誤差は明らかに独立ではありません．x, y がそれぞれ持つ独立な誤差に依存しています．このことより，楕円の当てはめに関しては最小二乗法と最尤推定法が異なります．
３．最尤推定法
じゃあ，最尤推定法はどうなるのか？という話ですが，工夫して式展開することによって見かけの確率変数を作り出します．まず，取得データの座標 (x と y の両方) に対して期待値０，分散が同一の独立な誤差が加わると考え，下記のようにあらわします．
上式を楕円の式に代入します．
上式を求める過程で，真値が楕円の方程式を完璧に満たすことと，誤差の二次の項は十分小さいとして無視しました．残った項を観察すると，楕円の方程式は下記の確率分布に従うことがわかります．
一点大事な点ですが，楕円の場合，分散が位置によって異なることがわかります．
上記の分散を共分散行列を用いて表現すると，下記のようになります．
確率変数を下記のように見なすと，
最尤推定法から求まる評価式 J は下記のように定義されます．この評価式のことをサンプソン誤差といいます．
よって，最尤推定法によって楕円当てはめを行う場合，上記評価式（サンプソン誤差）を最小にするようなパラメータセットを見つける必要がありますが，サンプソン誤差を最小にする方法としては，下記エントリで扱った Fundamental Numerical Scheme を用いて解くことができます．
daily-tech.hatenablog.com
Fundamental Numerical Sceme (FNS法)
ここでは，楕円の当てはめや基礎行列推定の際に評価式として出てくるサンプソン誤差を最小化する手段として，Fundamental Numerical Sceme を扱います．
まず，サンプソン誤差は下記のように定義されます．
この評価式は下に凸の関数であるので，最小値を持ちます．よって，パラメータ "θ" に関して微分し，微分値が 0 となるパラメータ を求めればこの値がサンプソン誤差を最小にする値になります．評価式の θ に関する微分は下記のようになります．
結局は上式，∇θJを０にするようなθを求めることに帰着します．但し，評価式は θ に関しての複雑な閉形式になっており，ストレートに ∇θJ = 0 を解くことができません．そのため，反復計算をします．また，上記の行列式が 0 になるということは，固有値が０となる固有ベクトルが存在するはずで，この固有ベクトルが求めるべき "θ" になります．アルゴリズムの手順を書き出すと．．．．
∇θJ = 0 の条件を満たすために，X を固有値分解して最も "影響の小さい" 固有ベクトルを取り出して X を再度固有値分解して，，，，というループでついには X の null space とそれを張る固有ベクトル (=解) が見つかるという流れですが，固有値分解を繰り返すと null space にたどり着くの下りがいまいちピンときてません．．．が，数学屋さんじゃないので次行きます．
最適化計算０ 最尤推定と最小二乗法
最適化計算の勉強をするにあたって，そもそも最小二乗法の下記の評価式
がどこから来るのかをまとめておきます．参考文献は金谷先生の教科書・PDFです．
２．尤度と最尤推定法
理論的には y = fθ(x) という数式にしたがう物理現象があるとします．ここで，関係式のパラメータセット θ が不明であり，実験を通してこれらのパラメータを取得しなければならなかったとします．この時，実験データに誤差が全くなければパラメータ個数分データを取得して方程式を解けばパラメータは求まります．しかし実際には実験データは誤差を含むので，複数のデータを取得して評価する必要があります．ここでは，複数データを用いて一番それっぽいパラメータを決定する方法として最尤推定法を導入します．
２．１ 誤差のモデル化
パラメータは数式を解いて求まるので，”誤差” をモデル化して数式に乗せてあげないといけません．そこで，下記のようにモデル化します．
ここで，誤差が期待値０，分散 σ2のガウス分布に従うとすると，誤差が ε になる確率は次のようになります．
実験データ一つ一つに対して理論式からのズレ（誤差）が定義できて，またそのズレ（誤差）が起こる確率は次のようになります．
２．３ 最尤推定法
N個の実験データを取得できているとすると，この”N個の実験データを発生させる確率が最も高いパラメータ”を求めるのが一番自然かと思います．この考え方を”最尤推定法”といいます．それぞれの誤差が生じる確率は数式としてわかっていますので，これを掛け算で求めていけば，N個の実験データが生じる確率が下記のように求まり，これを”尤度”といいます．
ということで，誤差が正規分布に従っていると仮定する場合，上式の
の部分を最小にしてあげれば確率が最大になることがわかります．着目点として，上記の評価式では分母の σ が定数であるため ∑ の外に出せます．このために，直線当てはめの場合はパラメータの正規化の方法によらず同一解が求まります．
３．最尤推定法を用いた最小二乗法（直線当てはめ）の導出
ここでは，２の最尤推定法の考え方を用いて直線当てはめの評価式を導出してみます．まず，理論式は下記のようになります．
ここで，数式 Ax + By + C の誤差は y と x 両方によって決まり，y と x のそれぞれの誤差を考慮すると，期待値０，分散 (A2 + B2) σ2の誤差が乗ることになります．つまり，最尤推定法から求まった式に代入すると下記のようになります．
ここで，直線当てはめに関しては最小化する評価式が上記のようになることがわかりました．直線当てはめの場合は，A2 + B2 + C2 = 1としても，A2 + B2 = 10 としても，A, B, C の定数倍に関しては不変であるので結果が変わらないことがわかります．つまり，下記二つの評価式のどちらを用いても，求まる結果は同じになります．
上記の説明で明らかになった通り，最小二乗法ではそれぞれの変数に対して誤差モデルに独立な正規分布を仮定し，最尤推定法を用いて取得データが起こる確率が最大になるようにパラメータを求めていることがわかります．
参考文献
参考文献１
- 作者: 金谷健一,菅谷保之,金澤靖
- 出版社/メーカー: 森北出版
- 発売日: 2016/10/27
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (1件) を見る
参考文献２
画像の三次元理解のための最適化計算［Ⅰ］
つくばチャレンジ２０１８参加チーム募集
ということで，下記のエントリ
daily-tech.hatenablog.com
を書いてちょうど一年．再びこの時期がやってきました．．．．そうです，つくばチャレンジ２０１８の参加募集開始です．（正確には 4/6 からエントリが始まるのでまだ参加募集は始まっていませんが．）
今年は３次元復元をマスターして，つくばチャレンジに生かすという大目標があります！が，昨年と同じくなかなか思うようには進まないもんですね．．．ただ，今年は４回目の参加ということで，完走します (ﾟДﾟ)ﾉ（と，言わせてください．）
今年は下記の３つをやります！
１．ハードを変える．
今年で４回目の参加になりますが，ハードを変えます！自分の作りこみが甘かったせいか，接触が悪くて止まったりということが頻発して時間を取られてしまうことが多かったので，思い切ります．で，さすがに今から自作は無理なのでつくばチャレンジ関係者の方が立ち上げてくれている下記サイトの i-Cart mini を使います．自作ロボットをやめるのは，，，，少し寂しいですが，しかたないですね！
t-frog.com
２．レーザースキャナを使う．
こちら， SI2017 で北陽電機の方とお話しさせていただく機会があり，なんとレーザスキャナを貸していただけることになりました．（北陽電機さん，本当にありがとうございました
www.hokuyo-aut.co.jp
３．画像ベースの３D地図を作成し，カメラで自己位置推定する．
あくまでも「画像処理だけで完走する！」という当初からのコンセプトは変えずこれからもやっていこうと思っていますが，今年は前年作った画像処理ベースの自己位置推定をもう少し発展させます．具体的にはコースの写真を一眼レフで撮影して３D地図を作成し，ロボットが自走するときにその地図のなかでどこを走っているかを判断させます．今やっている行列の勉強とかもこのためにやっている感じですが，ロボットに応用するところまで持っていけると今年の成果としては花丸です！
まずは試走会一回目（６／３０）の基本走行完了をマイルストンとして作業を進めていきます！
ということで，皆様応援よろしくお願いします！！！
最適化計算１ 直線の当てはめ 実験編
ということで，下記で理論をまとめたので直線当てはめの実験です．
daily-tech.hatenablog.com
ステップとしては７つあって，下記のステップに沿って実装してます．x
０．直線の係数定義
１．プロット用のデータを生成する．
２．フィッティングする誤差混入データを生成する．
３．生成したデータに対して，正規分布ベースの誤差をのせます．
４．最小二乗法を適用するマトリクス M を計算します．
５．固有値問題を解いて，パラメータベクトルを求めます．
６．最小二乗法で求まった直線（推定直線）を使ってデータを生成します．
７．推定直線のデータをプロットします．
ということで，結果は下記のようになりました．青が正しい直線，緑が推定直線です．当然ながら，誤差の分散を大きくしていくと青と緑の直線が合わなくなってきます．
※一点実装上の注意として！
２日くらいはまったのですが，Ax + By + C = 0 で求める係数のスケールが大きく違う (C >> A や C >> B, C << A や C << B) というような状況になると，丸め誤差の影響が顕著になってきてうまく推定直線が求まりません．これを回避するために，M行列を作成するときにスケーリングを行っています．
誤差分散 0.5
誤差分散 10
下記，実験に使った全コードです．
# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
def generateVecFromLine(a, b, c, minVal, maxVal, period):
nCols = (maxVal - minVal) / period + 1
if abs(a) > abs(b):
yVec = np.linspace(minVal, maxVal, nCols)
yVec = np.reshape(yVec, (yVec.shape[0], 1))
xVec = - (b / a) * yVec - (c / a)
else:
xVec = np.linspace(minVal, maxVal, nCols)
xVec = np.reshape(xVec, (xVec.shape[0], 1))
yVec = - (a / b) * xVec - (c / b)
return np.matrix(np.concatenate((xVec, yVec), axis=1))
def addGaussError(data, avg, stdv):
dataWError = data + np.random.normal(avg, stdv, data.shape)
return dataWError
def generateMMatrix(data):
xLen = max(data[:, 0])[0, 0] - min(data[:, 0])[0, 0]
yLen = max(data[:, 1])[0, 0] - min(data[:, 1])[0, 0]
factor = max(xLen, yLen)
a = np.matrix(factor * np.ones((data.shape[0], 1)))
P = np.concatenate((data, a), axis = 1)
return P.T * P, factor
def calcLineParamsFromEigenVector(M):
lamdas, v = np.linalg.eigh(M)
print(v[:, np.argmin(lamdas)])
return v[:, np.argmin(lamdas)]
if __name__ == "__main__":
print("Line Fitting Sample")
# 0. 直線の係数定義
A = -0.01
B = 0.025
C = 1
# 1. 直線上のデータを生成する．（真値）
dataOrg = generateVecFromLine(A, B, C, -100, 100, 0.01)
# 2. 誤差混入用データを生成する．（真値）
dataSampled = generateVecFromLine(A, B, C, -50, 50, 10.0)
fig, ax = plt.subplots(ncols = 1, figsize=(10, 10))
plt.xlim(-100, 100)
plt.ylim(-100, 100)
ax.plot(dataOrg[:, 0], dataOrg[:, 1])
# 3. 真値に対して，誤差を付与する．
dataWError = addGaussError(dataSampled, 0.0, 5.0)
ax.scatter(dataWError[:, 0], dataWError[:, 1], s = 2, edgecolors="red")
# 4. 行列 M を作成する．
M, factor = generateMMatrix(dataWError)
# 5. 固有値問題を解く
vEig = calcLineParamsFromEigenVector(M)
# 6. 最小二乗法で求まった直線からデータ生成
dataEstEig = generateVecFromLine(vEig[0, 0], vEig[1, 0], vEig[2, 0] * factor, -100, 100, 0.01)
# 7. 推定直線をプロット
ax.plot(dataEstEig[:, 0], dataEstEig[:, 1])
fig.show()
print("vEig : ")
print(vEig)
print("M : ")
print(M)
print("Factor : ")
print(factor)
最適化計算１ 直線の当てはめ 理論編
最近の投稿，行列やら最適化やら数学ネタ中心の投稿が続いてますが，これ全部「三次元復元」の理解のためです．おもった以上に数学のハードルが高く，一度きっちりと勉強しようと思いやり始めたらまあ時間がかかること．．．ここからは最適化の勉強をしていきます．今回の投稿では，ノイズを含むデータに対して直線を当てはめるという一番基本的なケースですが，主に金谷先生・菅谷先生の本・論文を拾い集めてきてまとめます．
１．直線当てはめ問題
実験で取得したデータをプロットしてまとめることって，高校とか大学の実験で結構やったと思います．実験なので所得したデータが何か既知の理論式に支配されているという前提で進めると思うのですが，取得データにはもちろん誤差が含まれるので完全にきれいには直線に乗りません．この時に，エクセルとか使って「一番それっぽい直線」を計算して求めたと思います．この「一番それっぽい直線」を求めるのが直線当てはめ問題です．
２．一番それっぽい直線とは？？
問題設定で「一番それっぽい直線」と書きましたが，「一番それっぽい」の定義を決めてあげないといけません．この種の直線当てはめ問題では，最小二乗法がよくつかわれますが，最小二乗法で使う評価式は下記になります．
上記の評価式ですが，実験で取得したデータが理論的には何かの直線上 "Ax + By + C = 0" に存在するとしたときに， それぞれの取得データ（点）と直線の距離が最も小さくなるような直線を求めています． ちなみに，直線を "y = ax + b" として最小二乗法を計算すると，評価式は下記のようになります．
この評価式では，評価データと直線の距離が最も小さくなるような直線ではなく，理論式と実験データの y の値のずれが最も小さくなる直線を計算しています．x を実験時に比較的正確に入力できるような場合はこちらのほうが感覚的には有用な気がしますね．
３．計算方法
計算方法ですが，自分が勉強した感じ，１．評価式の微分を 0 として計算する方法と２．評価式を二次形式に変換して固有値問題に帰着させる方法の二つがありましたが，ここでは２の線形代数アプローチで計算方法を求めてみます．
４．最小２乗法
つまり，，，，整理して書くとこの問題は下記のようになります．
参考文献
参考文献１
- 作者: 金谷健一,菅谷保之,金澤靖
- 出版社/メーカー: 森北出版
- 発売日: 2016/10/27
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (1件) を見る
参考文献２
画像の三次元理解のための最適化計算［Ⅰ］
