光学設計者の学習メモ: 機械学習

ラベル 機械学習 の投稿を表示しています。すべての投稿を表示

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回くらいで十分フィットしているように見えた。

2014年9月8日月曜日

Juliaで簡単な回帰

Juliaの練習がてら簡単な線形回帰を行った。

以下のことが分かる練習にはちょうど良いシンプルなコードだ。

乱数をどう発生させるか
計画行列をどう作るか
Gadflyを使って2つのプロットを重ねるにはどうするか
Gadflyを使って軸の設定をどうするか

これらがわかれば自作する当面の計算にはこまらない気がする…

linreg()はバイアス項を勝手に計算してくれるみたいなので、入れなくてよいみたいだ。

# #Regression Example

using Gadfly
using Distributions

N = 15
xdata = (rand(N) - 0.5) * 2 * pi
σ² = 0.5
function noisy(x, r=rand(Normal(0, σ²))
    return sin(x) + r
end

#X results in Nx9 array
t = noisy(xdata);
X = Float64[x^i for x=xdata, i=[1:9]] #bias is not needed for linreg(). it will interpolate automatically

#regression by linreg
coeffs = linreg(X, t)

#plotting
f = x -> coeffs[1] + sum([coeffs[i]*x^(i-1) for i = [2:10]])

x_reg = linspace(-pi, pi, 100)
y_reg = [f(xi) for xi=x_reg]

plot(layer(x=x_reg, y=y_reg, Geom.line), layer(x=xdata, y=t, Geom.point),
    Coord.Cartesian(xmin=-pi, xmax=pi, ymin=-1.5, ymax=1.5))

2014年8月10日日曜日

今さら聞けないRichardson-Lucy入門

デコンボリューションがとても魅力的だ。これを使えばレンズの収差やボケ、手ブレによる画像劣化が除去できるという夢の技術だ。

いろんな論文を読むとかならずといっていいほど以下の2つの手法が挙げられている。

1.Wiener Filter

Wikipediaに導出方法が載っている。

2.Richardson-Lucy Deconvolution

こちらもWikipediaに載っている。Wikipediaを読んでもわからなかったので、勉強した内容を少し解説してみる。

わかりやすいように1次元で話を進める。

をピクセルiの真の値(全くボケていないときの値)

をピクセルiで観測された値

とする。全くボケていないときはこの2つが一致する。でも今回はPSF(Point Spread Function)によってボケる。ボケがあるとピクセルjへ行く光は広がりをもってしまうので、ピクセルiへ少しおすそ分けしてしまう。
本来jへいくはずの光のボケでのih量を

とする。Richardson-Lucyではこのp(i,j)は既知だとするイメージは下の図。

ここでZ(i,j)なるものを導入する。Z(i,j)は期待値がP(i,j)λjのポアッソン分布に従いう確率変数。

これはピクセルiに届く光が確率的に変化し、期待値はp(i, j)λjということを意味する。

Wikipekiaにこの量は出てこないけど、元論文などを読むと出てくる。
もともとデコンボリューションはCTとかで研究されていて、CTは原子の崩壊(ガンマ線?)を捉えるので、観測される値が確率的になる、という経緯がある。
(一様な確率で起こる事象がある期間で観測される回数はポアッソン分布になる。詳しくはWikipediaをみてください。)

さて、点が1点だけなら復元も簡単かもしれんけど、周りの点の広がりが重なって観測されるので問題が難しくなる。例えば下の図のように重なるのでピクセルiで値は, 周辺のピクセルからの寄与を足しあわせたのものになる。

ここでzは確率的な値なので、観測結果と期待値から予測することで話を前に進めよう。
ピクセルiで観測された値をつかうと、その内ピクセルjからきた値は、割合で考えると下のような式になる。右辺の分母がiへくる光の総和で、分子がその内jからきたものを表す。

zはもともとjにあったものが周辺のピクセルiに散らばったものとかんがえることができるので、集めてやると真の値λjが分かる。真のzはわからないので推定値を使うとλを推定できてる。

この式にさっきのzの推定値を代入すると

このλjは真の値なので、知り得ない。そこで前回推定した値を使って逐次的に求めることにすれば

というふうになる。

実際に実装するときはWikipediaのようにコンボリューションの関数を使って書いたほうが楽だった。

(2014/11/15　式間違ってましたので修正しました。前に書いてたやつではPSFが対称じゃないとうまくいかなかったorz)

参考
https://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution
http://people.csail.mit.edu/dgreensp/may/deconvlucy.html

2014年6月20日金曜日

これなら分かる、最尤推定と事後確率最大化(MAP推定)の違い

最尤推定と事後確率最大化。どっちもガウス分布を仮定して対数を取ると2次式になるのでいままでなんとなく「同じようなもん」だと思っていた。

ここでは、観測された値をy、潜在変数をxとして話を進める。

まずはじめに訴えたいのは、どちらの手法も潜在変数xを求めたいということ。

尤度とは、あるモデル(仮説)を仮定したときの観測結果の確率P(y|x)である。xはいろいろな値が考えられるけど、観測結果を最もよく表すようにxを決定すれば現実に即していると言えよう。つまり、P(y|x)をxについて最大化するのが最尤推定。

次に事後確率最大化を考える。

まず単に事後確率といったらP(x|y)なのかP(y|x)なのかよくわからない。どちらも高校でならう事後確率の形式やし。

ここでは、観測がされた後の確率を事後確率と呼ぶ。観測が条件となる確率。つまりP(x|y)のこと。P(x|y)を最大化するの事後確率(MAP推定)だ。では最尤推定とどうちがうのか??

ベイズの定理を使うと

P(x|y) ∝P(y|x)P(x)

右辺は(尤度)*(事前確率)となっている。

ここで(事前確率)を何かの定数だと考えてP(x|y)を最大化を考えると、上の尤度を最大化する場合と同じになっている!

つまり、事後確率最大化は最尤推定に事前確率というバイアスがかかったものだということだ。事後確率は主観とか言われてこれがベイズと呼ばれるアプローチだ。

事後確率の負の対数を取って最大化すると正則化項として出てくる(よく教科書にでてくるパターン)。つまり主観がある種のブレーキとして働くと考えられる。

2014年6月3日火曜日

scikit-learnを使って数字認識(2)　k-NNを使った学習

scikit-learnを使えば簡単な機械学習ならすぐにできる。
今回は特徴ベクトルとして画素値をもちい、それをk-NearestNeighborsにぶちこむ。
k-meansもためしたけど、こちらはあまりうまく行かなかった。
Pythonだと学習した識別器(オブジェクト)をPickleで保存して簡単に取り出せる。

10種類くらいのフォントを学習させて11種類目でためしたが、正答率は100%だった。今後もっといろんな状況下(画素が少ないサンプルとか)で試していってどこまでいけるか確かめたい。

以下メモ

学習部分

import numpy as np
import cv2
from sklearn.neighbors import KNeighborsClassifier
from sampling import convert_to_binary
"""
1. Read sample image and convert to 1d feature vector
2. pass the feature vectors to kmeans clustering maching
3. pickle the result.
"""

def convert_to_feature_vector(src):
    """
    convert source image to a 1d feature vector
    """
    N =10
    im = cv2.GaussianBlur(src,(5,5),0)
    im = cv2.resize(src, (N, N))
    im = im.reshape(N**2)
    #im = np.array(im>124, dtype=np.int8) #convert to 0 and 1
    return im


if __name__ == "__main__":
    #read jpgs, and resize them as a vector

    dirname = [str(i) for i in range(10)]
    dirname += ["dot"]
    dirname += ["bar"]

    X = [] #sample data
    Y = [] #label
    for dn in dirname:

        #input label(0-11)
        if dn == "bar":
            label = 11
        elif dn == "dot":
            label = 10
        else:
            label = int(dn)


        fnames = os.listdir(dn)
        fnames.remove(".DS_Store")

        for fn in fnames:
            #print os.path.join(dn, fn)
            im = cv2.imread(os.path.join(dn, fn))
            im = cv2.cvtColor(im, cv2.COLOR_BGR2GRAY)
            im = convert_to_feature_vector(im)
            #print im.reshape(10,10)
            X.append(im)
            Y.append(label)

        #convert to np.array
    X = np.array(X)

    #
    knn = KNeighborsClassifier(n_neighbors=5)
    knn.fit(X, Y)

    #save the classifier as pickle
    import pickle
    with open("knn_trained.dump", "w") as f:
        pickle.dump(knn, f)

予測部分

#coding: utf-8

from sklearn import cluster
import pickle
import cv2
from training_knn import convert_to_feature_vector


if __name__ == "__main__":
    with open("knn_trained.dump") as f:
        knn = pickle.load(f)
        print knn

    import os
    dirname = "test_data"
    fnames = os.listdir(dirname)
    try:
        fnames.remove(".DS_Store")
    except ValueError as e:
        print e


    for fn in fnames:
        im = cv2.imread(os.path.join(dirname, fn))
        im = cv2.cvtColor(im, cv2.COLOR_BGR2GRAY)
        x = convert_to_feature_vector(im)
        y = knn.predict([x])
        print fn, y

scikit-learnを使って数字認識(1)　サンプル取得

OCRはフリーのものもあるけど、数値データを起こしたいときには無駄に高機能になってしう。それだけでなくご認識も増える。

手書きでない数値(0-9, "-", ".")を認識するだけならそんなに苦労しないのでは?と思ったのでちょっと実験してみた。

フツウの人は手書き文字認識に興味があるやろけど、
調べてみるとStackoverflowに以下のような投稿があって結構盛り上がっている。
https://stackoverflow.com/questions/9413216/simple-digit-recognition-ocr-in-opencv-python

今回は

上のリンクの方法で数値を切り取り
scikit-learnを使って学習
学習データを元に新しい予測

ということを行う。

まず、以下のように切り取られた画像から数値を取得したい。

イメージしている成果物は、WebやPDF上をマウスで矩形選択するとそこの数値をテキストデータに変換するツールだ。

ほぼ上のサイトと同じことをしているが、一つだけ違うところfindContoursで輪郭の取得オプションだ。opencvでは輪郭の取得のみならず、それら輪郭の関係まで返してくれる。今回はcv2.RETR_CCOMPを指定した。(一番外側の輪郭とするとなぜか図の枠が検出されたので)輪郭についてはこのページがわかりやすかった。輪郭の階層を指定することで輪郭の大きさによる分類など無駄な作業が省ける。

下のプログラムを実行するとこうなる。

if __name__ ==  "__main__":
    sample_dir = "fonts"
    sample_file = "sample03.png"
    im = cv2.imread(os.path.join(sample_dir, sample_file))
    im_copy = im.copy()
    gray = cv2.cvtColor(im,cv2.COLOR_BGR2GRAY)
    blur = cv2.GaussianBlur(gray,(5,5),0)
    thresh = cv2.adaptiveThreshold( src=blur, 
                                    maxValue=255, 
                                    adaptiveMethod=cv2.ADAPTIVE_THRESH_MEAN_C,#cv2.ADAPTIVE_THRESH_GAUSSIAN_C,
                                    thresholdType=cv2.THRESH_BINARY,
                                    blockSize=5,
                                    C=3)
    thresh_copy = thresh.copy()  #thresh are destroyed when findCountours

    cv2.imshow('threshold', thresh)
    cv2.waitKey(0)


    contours,hierarchy = cv2.findContours(thresh_copy, cv2.RETR_CCOMP, cv2.CHAIN_APPROX_SIMPLE)
    print "%d contours are found" % len(contours)


    i = 0
    for con, hie in zip(contours, hierarchy[0]):
        if hie[3] != -1: #2:first child, 3:parent
            (x,y,w,h) = cv2.boundingRect(con)
            #roi = thresh[y:y+h,x:x+w]
            roi = im[y:y+h, x:x+w]

            cv2.rectangle(im_copy,(x,y),(x+w,y+h),(0,0,255),2)
            cv2.imshow('rect',im_copy)
            #cv2.imshow("threshold", roi)
            
            key = cv2.waitKey(0)
            if key == 27: # (escape to quit)
                sys.exit()

            c = chr(key)

            if c == ".":
                c = "dot"
            elif c == "-":
                c = "bar"
            elif c == "n": #negative
                c = "negative"

            file_head = sample_file.replace(".png", "")
            if c == "negative":
                fname = os.path.join(c, file_head + "_" + c + ("%02d.png"%i))
                i+=1
            else:
                fname = os.path.join(c, file_head + "_" + c + ".png")

            #write ROI 
            print fname
            #cv2.imwrite(fname, roi)

登録: 投稿 (Atom)