syuntoku14の進捗

進捗を書きなぐります

進捗(5/24) 病み上がり

 

進捗(5/24) 確率的プログラミング楽しいよね

病み上がりなので特に何もしていない。

ipad環境を構築したりしていた。ipad非常に良いですね。あとはペーパーライク画面保護フィルムがあれば完璧。

ブックスタンドも導入した。クソ便利になった。学習環境がとても充実している。

 

Coursera

Why model p(x)

  • Generate new data
  • Detect anomalies and outliers
  • Work with missing data
  • Represent your data in a nice way

How to model p(x)

  • \(\log \hat p(x)=CNN(x)\)

However, in this case we have to compute \(p(x)=\frac{\exp (CNN(x))}{Z}\)

which is expensive to compute when x is really high dimention. So this is infeasible.

  • Use the chain rule

\[p(x_1,...,x_d)=p(x_1)p(x_2\mid x_1)...p(x_d \mid x_1, ....., x_{d-1})\\ =RNN(x_1,...,x_{k-1})\]

Reference: "Pixel recurrent neural networks"(2016)

It is more easier to compute because now we have to consider the normalization constant which is based on one dimentional distribution.

However, this is very slow.

  • \(p(x_1,...,x_d)=p(x_1)...p(x_d)\)

But this assumption is too restrictive. Each pixel must have something to do with each other.

  • Mixture of several Gaussians(GMM)

Still too restrictive

  • Mixture of infinitely many Gaussians

\[p(x)=\int p(x\mit t)p(t)dt\]

The idea is that the image x has a corresponding latent variable t, and the image x is caused by the t. The conditianl distribution of \(p(x \mid t)\)is Gaussian.

Using CNNs with a mixture of Gaussians

  • Continuous mixture of Gaussians

\[p(x)=\int p(x\mid t)p(t)dt\]

\[p(t)=N(0,I)\]

\[p(x\mid t)=N(\mu (t), \Sigma(t))\]

if $ (t) = Wt+b, (t) =_0$ get PPCA

But if x is image, why not \[\mu(t) = CNN_1(t)\]

\[\Sigma(t) = CNN_2(t)\]

However, if our image is 100x100, our covariance matrixes is 10000 x 10000. This is really large.

To avoid this, let's say our covariance will be diagonal.

Learning

\[\max_w p(X\mid w)=\int p(X\mid T,w)p(T)dt\]

This is Latent Variable model... But we can't use normal EM algorithm.

Need to compute \(p(T\mid X,w)\)

\[\log p(X\mid w)\geq L(w,q)\]

\[maximize_{w,q} L(w,q)\]

In this step, we have to compute the integral and the integral conains CNN so this hard to do analitically.

How about MCMC?

\[\mathbb{E}_q \log p(X,T\mid w)\approx \frac{1}{M}\sum^M_{s=1}\log p(X,T_s\mid w)\\ T_s\sim q(T)\]

In this method, we have to wait the MCMC conversion every time of update. This is very slow.

We can try variational inference.

\[log p(X\mid w)\geq L(w,q)\]

\[maximize_{w,q} L(x,q)\]

subject to \(q_i(t_i)=\hat q(t_{i1})...\hat q(t_{im})\)

But again intractable.

 

pythonで体験するベイズ推論

教科書p.15 PyMCの入門

例として、日々のメール受信数をモデリングすることを考えよう。i日目のメッセージ数を\(C_i\)とし、これがポアソン分布に従うと仮定する。

\[C_i \sim Poi(\lambda)\]

ここで、ポアソン分布は

\[P(Z=k)=\frac{\lambda^ke^{-\lambda}}{k!}\]

で、確率変数Zがポアソン分布に従うことを \[Z \sim Poi(\lambda)\]

と書く

話に戻り、メッセージ数のモデルを考えよう。 ここで、メッセージ数の急激な変化点において、\(\lambda\)が変化していると仮定する。つまり、 \[\lambda = \left\{ \begin{array}{ll} \lambda_1 \;(t<\tauのとき)\\ \lambda_2 \;(t\geq \tau のとき) \end{array} \right. \]

ここで、\(\lambda\)自体が指数分布に従っているとする。

つまり

\[\lambda \sim Exp(\alpha)\]

ここで、\(\alpha\)を係数データの平均を取って、その逆数を取るという手法で定義しよう。

\[\frac{1}{N}\sum^N_{i=0}C_i\approx E[\lambda\mid \alpha]=\frac{1}{\alpha}\]

\(\tau\)の決め方を考えよう。データには結構なばらつきがあるので、一様分布を用いるのがいいだろう。

\[\tau \sim DiscreteUniform(1,70) \\ \rightarrow P(\tau = k)=\frac{1}{70} \]

この未知変数の事前分布は、結構複雑な形だ。ほしいのは事後分布であって、あまり事前分布に気を向けたくない。そこで、PyMCを使ってこれを解いてみよう。

モデルは次のようにコードに起こすことができる

import pymc3 as pm

with pm.Model() as model:
    alpha = 1.0 / count_data.mean()
    
    lambda_1 = pm.Exponential("lambda_1", alpha)
    lambda_2 = pm.Exponential("lambda_2", alpha)
    
    tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data-1)

これらはPyMC3のstochastic variablesと呼ばれる。stochastic変数はバックエンドで乱数生成器として扱われるのでそのような名前がついており、以下のようなrandomメソッドで呼び出すことができる。

print("Raondom output:", tau.random(), tau.random(), tau.random())

次に、lambda_1とlambda_2をスイッチングさせる関数lambda_を考えよう。これは、実際には確率変数と考えることができる。

with model:
    idx = np.arrange(n_count_data)
    lambda_ = np.math.switch(tau > idx, lambda_1, lambda_2)
with model:
    observation = pm.Poisson("obs", lambda_, observed=count_data)

ここでは、データ生成モデルであるポアソン分布のオブジェクトを生成し、それを変数observationが受け取っている。observed=count_dataにすることで、値が観測値であり、解析時には固定されていることをPyMCに伝えている。

with model:
    step = pm.Metropolis()
    trace = pm.sample(10000, tune=5000, step=step)

上のコードは学習ステップとみなすことができる。メトロポリスヘイスティング法によるMCMCを行い、事後分布を求めている。

求まったやつ:

ここからわかること: * \(\lambda_1\)はおおよそ18で、\(\lambda_2\)はおよそ23である。二つの分布は明らかに異なっているので、ユーザのメッセージ受信に明らかな変化があったことがわかる。

  • 事後分布が明らかに指数分布ではない。この事後分布はここに登場するどの分布とも異なっている。 数学的な解析を行った場合、数式が複雑すぎて手に負えなくなってしまっていただろう。
  • \(\tau\)の分布を見ると、離散確率変数なので、事後分布が確率質量関数である。45日めにユーザが振る舞いを変えた可能性が50%である。

最後に事後分布からのサンプルを使って、t日目のメッセージ数の期待値はいくらか?といおう問題に答えてみよう。ポアソン分布の期待値はパラメータ\(\lambda\)なので、t日目における\(\lambda\)の期待値は何か?という問題と同じになる。

# tau_samples, lambda_1_samples, lambda_2_samples contain
# N samples from the corresponding posterior distribution
N = tau_samples.shape[0]
expected_texts_per_day = np.zeros(n_count_data)
for day in range(0, n_count_data):
    ix = day < tau_samples  # ix is like [True, False, True, True...]
    expected_texts_per_day[day] = (lambda_1_samples[ix].sum()
                                   + lambda_2_samples[~ix].sum()) / N

上のコードで、lambda_1_samples[ix].sum()/Nは\(\lambda_1\)のt日目の期待値を表している。

同様に\(\lambda_2\)についての期待値も計算して和を取れば、t日目のメッセージ数の期待値が求まる。

 

ROSで自作ロボットのシミュレーション

ROSでは、UDRFという形式で物体のハードを記述する。これがとにかくめんどくさい。やりたくない。