光学設計者の学習メモ: 9月 2015

2015年9月23日水曜日

ipython notebookでVim(バインディング)

ipython notebook(Jupyter)でもvim bindingが使えるらしい。
https://github.com/lambdalisue/jupyter-vim-binding

ipython notebookで使いにくかった点が1つ改善された。

環境はMacにAnacondaを導入している。
基本的な導入は上のGithubのリンクに書いてあるけど…一応メモしておこう。

まずはipythonをアップデートして、さらにjupyterを入れる。
(Jupyterを入れないとnotebookモジュールがないと怒られる)

$ conda update ipython
$ conda install jupyter

次に以下を実行

from notebook.nbextensions import install_nbextension
install_nbextension('https://goo.gl/5TK96v', user=True, destination="vim_binding.js")

一時的に使うときは

%%javascript
Jupyter.utils.load_extensions('vim_binding')

ずっと有効にしたいときは

%%javascript
Jupyter.notebook.config.update({
  'load_extensions': { 'vim_binding': true },
});

これでipython notebook上でVimバインディングが使える。いやー、便利になった。

2015年9月13日日曜日

Iteratively Reweighted Least Squares　についてサクッと。

Iteratively Reweighted Least Squares(IRLS)ってたまに出てくるけど、何か分かってなかった.

Least Squares

LSでは線形モデル、

$\bf{y} = \bf{\beta}^T \bf{x} + \bf{\epsilon}$

というモデルを考えたときにデータとモデルの2乗和誤差を最小にするように $\bf{\beta}$ を決定する。このとき誤差 $\bf{\epsilon}$ を1つの定数で書ける(データのばらつきが次元で同じ)ときは擬似逆行列を使って係数 $\bf{w}$ を決定することができる。導出はないけど、微分が0になるように $bf{\beta}$ を決定すれば

$\bf{\beta}^{*}= (X^TX)^{-1}X^T \bf{y}$

Weighted Least Squares

先程はデータの次元でばらつきが同じとしたけど、次元間の分散が違うことだってある。そんな場合にWLSをつかう。調べた感じWikipediaが一番分かりやすい。
https://en.wikipedia.org/wiki/Linear_least_squares_(mathematics)#Weighted_linear_least_squares
この場合、

$\bf{\beta}^* = (X^TWX)^{-1} X^TWy$

となる。ただし $W$ は対角成分に分散の逆数を並べた行列。例えば、あるデータ次元についてはほぼ確定的なら始めから分散は小さくできるので、それを反映させることができる。式を見ると、データを標準偏差で割っているのと同じになっているので、データが事前にある場合は標準化(standardization)しとけば同じことになりそうだ。

Iteratively Reweighted Least Squares(IRLS)

上の2つは一発で最適な係数 $\bf{\beta}^*$ を求めることができた。これは2乗和誤差が $\bf{\beta}$ の2次で書けるからだ。IRLSでは、2次でない場合でも2次で近似しながら重みWを更新して、最適な係数 $\bf{\beta}$ を求めてくれる。これまたWikipediaがわかりやすかった。
https://en.wikipedia.org/wiki/Iteratively_reweighted_least_squares

Lpノルムを最小化したいときは、
$\bf{\beta}^{(t+1)} = (X^TW^{(t)}X)^{-1} X^TW^{(t)}y$

$W^{(t)}$ は対角行列。 $W^{(0)}=\bf{1}$ と初期化しておき、対角成分の各要素は以下のように更新していく。
$w_i^{(t)} = | \bf{y}_i - X_i \beta^{(t)}|^{p-2}$

IRLSはスパースな係数を選択するためにL1最適化を使う場合に用いられたりするようだ。
コンピュータビジョン最先端ガイド6巻の3章に少しだけIRLSの記述がある。
Scikit-learnなどではOMPが実装されているみたいだが。

2015年9月1日火曜日

Python でPRML7.2章　RVMを実装。

いつ使うかは分からんが、RVMを実装しておこう。
7.2章の数式をそのまま書いた。動作確認は1次元のガウス基底で行った。

def rvm(t, PHI, alphas0, beta0, maxloop=10):
    """
    t: observed values.
    X: designe matrix
    alphas0: initial weights for each parameter.
    beta0: an initial precision value.
    
    returns
    alpha: prior weights
    beta: prior precision
    m: prior means
    V: prior covariance matrix
    
    See PRML 7.2 for details.
    
    ** to many loop causes 0 devision for some reason.
    """
    beta = beta0
    alphas = alphas0
    
    for i in range(maxloop):
        A = np.diag(alphas)
        V = np.linalg.inv(A+beta*np.dot(PHI.T, PHI))
        m = beta*np.dot(V, np.dot(PHI.T, t))
        gammas = 1 - alphas * np.diag(V)
        
        #new alphas and beta
        alphas = gammas / m**2 
        beta =  (N - np.sum(gammas)) / np.linalg.norm(t - np.dot(PHI, m))
        
    return alphas, beta, m, V

10個のガウス基底でフィッティングをした。
今回書いたコードではループの回数が多すぎるとzero divisionでエラーとなった。でもループは10回くらいで十分フィットしているように見えた。