2015-08-02

Marsaglia法とBox-Muller法

一様分布からガウス分布を作るアルゴリズムとして有名なBox-Muller法というのがあります．
式としては、 ${z_1,z_2\in(0,1)}$ とした時に
${\displaystyle y_1 = \sqrt{\mathstrut -2\ln z_1}\cos 2 \pi z_2 \\ y_2 = \sqrt{\mathstrut -2\ln z_1}\sin 2 \pi z_2 }$
という変換を行う．これがBox-Muller
ボックス＝ミュラー法 - Wikipedia

それと似たもので、一様分布する変数の範囲がちょっと変わって ${z_1,z_2\in(-1,1)}$ としてやって、更に
${ z_1^2+z_2^2 \leq 1 }$
を満たすものだけを取り出して（要するにそうじゃないやつを棄却して）
${ y_1 = z_1(-2 \ln r^2 / r^2)^{1/2}\\ y_2 = z_2(-2 \ln r^2 / r^2)^{1/2} }$
という変換をしてやるのが、Marsaglia
Marsaglia polar method - Wikipedia, the free encyclopedia

なんか読むとBox-Mullerのほうがcosとsinの計算を並列でできる？から早いとかなんとか．まあ棄却する部分があるMarsagliaのほうが遅そうやなぁという雰囲気はあるよね．

というわけで実装

import numpy as np
import math as math
import time
import matplotlib.pyplot as plt

pears = 1000
##compare marsaglia with box-muller
def boxmuller():
    """
    using box-muller method
    return:points random sampling following N(0,1)
    """
    points = []
    for i in range(pears):
        z1 = np.random.uniform(0,1)
        z2 = np.random.uniform(0,1)
        y1 = np.power((-2 * math.log(z1)),0.5) * math.cos(2 * np.pi * z2)
        y2 = np.power((-2 * math.log(z1)),0.5) * math.sin(2 * np.pi * z2)
        points.append([y1,y2])
    return points

def marsaglia():
    """
    marglia method
    return:points random sampling following N(0,1)
    """
    points = []
    for i in range(pears):
        z1,z2 = 0.,0.
        r = 0.
        while 1:
            z1 = np.random.uniform(-1,1)
            z2 = np.random.uniform(-1,1)
            r = z1 * z1 + z2 * z2
            if r <= 1:
                break
        y1 = z1 * np.power(-2. * math.log(r) / r,0.5)
        y2 = z2 * np.power(-2. * math.log(r) / r,0.5)
        points.append([y1,y2])
    return points

methods = [boxmuller,marsaglia]
for i in methods:
    start_time = time.time()
    for k in range(10):
        i()
    end_time = time.time()
    print (end_time-start_time)/100,i

points = np.array(marsaglia())
plt.plot(points[:,0],points[:,1],"ro",alpha=0.5)
points = np.array(boxmuller())
plt.plot(points[:,0],points[:,1],"go",alpha=0.5)
plt.show()

結果

0.00240999937057 <function boxmuller at 0x04357D70>
0.00396000146866 <function marsaglia at 0x04357DB0>

やっぱboxmullerのほうが早いみたい．

f:id:dette:20150802093017p:plain

どっちも標準正規分布にはなってるみたい．

2015-07-28

PRML 第5章

ニューラルネットワークのところを読んでいますが詰まりました．

ニューラルネットワークでは、最終的な出力が

${\displaystyle y_k(x,w)=\sigma(\sum^M_jw^{2}_{ik}h(\sum^D_iw^{1}_{ji}x_i)) }$

という風になっていて、

そこで疑問なのが、
例えばKクラスの分類問題をニューラルネットワークに説かせようと思った時、
この最後の ${\sigma}$ の部分をソフトマックス関数

${\displaystyle \sigma=\exp(a_k)/\sum_j\exp(a_j) }$
とすれば良い
（ただし ${a_k}$ は出力ユニットへの活性である）

という説明があるのだけれど、疑問がひとつ．

最終の出力ユニットの活性 ${a_k}$ は ${a_k=\sum^M_kw^{2}_{kj}a_j}$ となっているのだから、
出力ユニットjへの入力は既に隠れユニットの出力 ${a_j}$ と重み ${w_{kj}}$ の積の和になっているはずで、

そこからどうやって各々の ${a_j}$ を計算しなおしてるんだろうか…という

なんか多分どうでもいい勘違いなんだろうけど凄くイライラ．

2015-07-27

prml 第4章

PRML

先週から読み始めてようやく第四章までたどり着きました。噂には聞いていたのですが、ここまで重たいものだとは正直予想外で、なんとかかんとか頑張ってついていっているのか、はたまたわかったような気になっているのかは謎ですが、とりあえずちょっとづつ前進している感じです。

読んでいて思うのは、式変形はわかるんだけど、その意味合いを考えるとちょっとこんがらがるというか、冷静に考えると意味わかってないなって気づくところが多いということ。これはベイズの定理を使うのに慣れ始めてから少しづつ改善はしているのですが、やっぱりなぜこの式をここで使うのか？とかをちゃんと考えながら読まないと、ビショップさんの意図を汲み取れなくなっていくので、そこがしんどいです。

あと自分が理解できていないからかも知れませんが、微妙に天下り的に結果を先に与えていたり、すでにわかっているパラメータを条件付き確率の既知の部分から省いていたり、そういう微妙な罠がそこかしこにあって、記述は綺麗だけど意味を汲み取るのに時間が掛かる時があって、そこも苦労している点です。

まあわかってる人からすれば大したことはない省略だとは思うんですけども。モデルの選択の部分でモデル選択のMiが突然なくなった時とか、結構悩まされました。みんなどうやってクリアしてるのか気になるところです。

やっぱ誰か一緒に読む人がいたほうがいいなぁ…

あと、今まで研究室にあったやつを借りて読んでいたのですが、やっぱ高いけどこの本は買わないとダメだと決心がついたので、そろそろ買いに行こうかと思います。買ったら投げ出さんやろうという、退路を断つ作戦。そうでもしないと最後まで絶対やらんしな。

2015-07-23

カーネル密度推定

PRML 機械学習 Python

PRMLを読んでいます。とりあえず出てくる式はすべて自分で追って、変形とかもやってるのでなかなか進まなくてもどかしいですが、もう多分確率密度関数とかを復習することは無いだろうし…と思いながらやってます。

今日は二章のおわりまで行きました。
ノンパラメトリックの推定方法のところで、この間「外野フライをカーネル密度推定を用いて解析する」っていう論文を読んで気になってたカーネル密度推定のお話が出てきました。ざっくり言うと、その点の周りにいっぱい点があれば、そこの確率は高いよねっていうことです。発想としては至ってシンプルなんですね。直感的。

なんか可視化してみたら面白そうだったので、pythonの練習も兼ねて実装をしてみました。
やってるのは、一次元のガウス分布に従うサンプルから確率密度関数を推定するっていうものです。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal

class KernelEstimation():
    def __init__(self,trainingdata,h,D,kernel):
        self.trainingdata = trainingdata
        self.h = h
        self.D = D
        self.kernel = kernel

    def labeltext(self):
        return "kernel: "+self.kernel+" sigma: "+str(self.h)

    def _kernel(self,d):
        """
        return kernel distance
        """
        if self.kernel == "gausian":
            norm = np.linalg.norm(d)
            return np.exp(-norm*norm/(2.*self.h*self.h))
        return 0
    def _normalize(self):
        """
        normalization parametor
        """
        if self.kernel == "gausian":
            return 1/np.power((2*self.h*self.h*np.pi),(self.D/2.))/self.trainingdata.size
        else:
            return 0


    def predict(self,x):
        """
        predict by kernal estimate
        """
        prob = np.zeros(x.size)
        for i in range(x.size):
            for k in range(self.trainingdata.size):
                prob[i] += self._kernel(x[i]-self.trainingdata[k])
        return prob*self._normalize()
    
rv = multivariate_normal(1,1)
x = np.linspace(-5,5,100)
plt.plot(x,rv.pdf(x),"k-")
samples = rv.rvs(100)#sampleの生成
plt.hist(samples,normed=True, facecolor='green',alpha=0.3)

kernel = KernelEstimation(samples,0.30,1,"gausian")
plt.plot(x,kernel.predict(x),"b--",label=kernel.labeltext())
plt.legend()
plt.show()

KernelEstimationというクラス作ってみました、完全に練習ですが少しは見通しが良くなってるかな？
f:id:dette:20150723191716p:plain
sigmaの値を変えると推定の雰囲気が色々と変わるので面白いですね。

nykergoto’s blog

機械学習とpythonをメインに

Marsaglia法とBox-Muller法

PRML 第5章

prml 第4章

カーネル密度推定