PyMCのお勉強

目次 (2023/4/11 ver.0.01)

1番目のタブが選択された！

1.はじめに

もくじ

1.1 NumPyroとは

1.NumPyroとは
  • PyMCとは簡単に言えばPython向けのベイズ推定用のサンプリングモジュールです。
  • 確率的プログラミング言語(PPL)とも呼ばれます。Python向けに限らず同等なものはいくつか存在しています。
  • Pyro, NumPyro, TensorFlowProbability, Stan(PyStan), WinBUGS, JAGSあたりでしょう。
  • どれもMCMCサンプリングと変分ベイズ(変分推論、ADVI)、プラスアルファで確率的最適化ができるはずです。
  • PyMCの特徴は開発が盛んなこと、そして使用者が多いことです。
  • 以下はGoogle Trendsで調べた上記のPython向けPPLの人気です。
2.データセットの作成
  • 今回は、scikit-learn の tutorial でも使用されている 以下の関数に対し GP モデルの構築を行います。
  • モデル構築用データは、x=(1,3,5,6,7,8) の6点です。

    • $\Large \displaystyle y=x\sin(x)$

以上

2.ガウス過程

もくじ

3.1 ガウス過程の例

本節は、「 ケモインフォマティクスのオンライン入門書 」を引用している。

1.データセットの作成
  • 今回は、scikit-learn の tutorial でも使用されている以下の関数に対し GP モデルの構築を行います。
  • モデル構築用データは、x=(1,3,5,6,7,8) の6点です。
2.GPモデルの構築
  • 目的変数の scaling を行った後に、GP モデルの定義および学習を行います。
  • 今回は kernel 関数に以下の構成を用いました。

    • $\Large \displaystyle (Constant Kernel)×(RBF Kernel)+(White Kernel)$

  • Constant Kernel が定数 kernel、RBF Kernel がRadius Based Function kernel (あるいは Gaussian kernel)を表します。
  • White Kernelは目的変数のノイズの大きさに関する kernel 関数で、(今回は関係ありませんが)データセットにノイズが含まれる場合に重宝します。
  • White Kernel を導入する場合、GaussianProcessRegressor の引数 alpha を0とすることに注意してください。
  • GaussianProcessRegressor の fit 関数を実行すると、GP モデルの学習が行われます。
  • 特に指定がなければ、モデル構築用データに対する尤度が最大となるように kernel 関数のパラメータが最適化されます。
3.GPモデルを用いた予測
  • 構築した GP モデルを用いて、各 x に対する y の値の95%信頼区間を予測します。
  • GaussianProcessRegressor の predict 関数は、return_std=True とすることで予測される目的変数の分布の平均および標準偏差を出力します。 
    • # coding: UTF-8
import numpy as np
from matplotlib import pyplot as plt
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import ConstantKernel, RBF, WhiteKernel
from sklearn.preprocessing import StandardScaler

# 解析関数の定義
def f(X):
    return X * np.sin(X)

# 1.データセットの作成
X = np.atleast_2d([1., 3., 5., 6., 7., 8.]).T
y = f(X)


# データセットの可視化
plot_X = np.atleast_2d(np.linspace(0, 10, 1000)).T
fig = plt.figure(figsize=(8, 6))
plt.plot(plot_X, f(plot_X), 'k')
plt.plot(X, y, 'r.', markersize=16)
plt.xlabel('$x$', fontsize=16)
plt.ylabel('$y$', fontsize=16)
plt.ylim(-8, 12)
plt.legend(['$y = x*\sin(x)$', 'Observed values'], loc='upper left', fontsize=16)
plt.tick_params(labelsize=16)
plt.show()

# データのscaling
# scikit-learnに実装されているStandardScalerを利用
# 説明変数のscalingはしなくても問題ありませんが、目的変数のscalingは必須です(平均の事前分布を与えるのが難しいため)
scaler_y = StandardScaler().fit(y)

# 2.GPモデルの構築
kernel = ConstantKernel() * RBF() + WhiteKernel()
gpr = GaussianProcessRegressor(kernel=kernel, alpha=0)
gpr.fit(X, scaler_y.transform(y))

# kernel関数のパラメータの確認
print(gpr.kernel_)
# 1.36**2 * RBF(length_scale=3.07) + WhiteKernel(noise_level=0.061)

# 3.GPモデルを用いた予測
# plot_Xに対する平均および標準偏差の予測
pred_mu, pred_sigma = gpr.predict(plot_X, return_std=True)
pred_mu = scaler_y.inverse_transform(pred_mu)
pred_sigma = pred_sigma.reshape(-1, 1) * scaler_y.scale_

# 各xに対する95%信頼区間を表示
fig = plt.figure(figsize=(8, 6))
plt.plot(plot_X, f(plot_X), 'k')
plt.plot(X, y, 'r.', markersize=16)
plt.plot(plot_X, pred_mu, 'b')
# データが正規分布に従う場合、95%信頼区間は平均から標準偏差の1.96倍以内の区間となる
plt.fill_between(plot_X.squeeze(), (pred_mu - 1.9600 * pred_sigma).squeeze(), (pred_mu + 1.9600 * pred_sigma).squeeze())
plt.xlabel('$x$', fontsize=16)
plt.ylabel('$y$', fontsize=16)
plt.ylim(-8, 12)
plt.legend(['$y = x*\sin(x)$', 'Observed values', 'Predicted mean', '95% confidence interval'], loc='upper left', fontsize=16)
plt.tick_params(labelsize=16)
plt.show()
  • 出力されたグラフを見れば分かる通り、モデル構築用データが存在する領域では95%信頼区間が狭く(予測された標準偏差が小さく)、モデル構築用データが存在しない領域では95%信頼区間が広く(予測された標準偏差が大きく)なっています。
4.GPモデルを用いた実験計画法
  • GP を用いて、簡単な実験計画法(Experimental Design)を行ってみましょう。
  • GPを用いた実験計画法は一般的にベイズ最適化と呼ばれ、少ない実験回数での設計変数最適化 [5] や機械学習モデルのハイパーパラメータ最適化 [6] に活用されています。
  • ベイズ最適化では、GP モデルで予測された目的変数 y の分布を元に、Acquisition Function と呼ばれる指標を計算し、Acquisition Function が最大となる点を次の実験対象として選択します [2]。
  • 目的変数の最小化を目指す場合、予測された y の分布の平均が小さく(予測値が小さい)、標準偏差が大きい(不確かさが大きい)ほど Acquisition Function は大きな値を取ります。
  • 構築した GP モデルを用いて、より小さい観測値が得られる説明変数 x の値を探索します。
  • 今回は Acquisition Function として、Expected Improvement (EI) を用います。
  • EI は以下の式で計算され、モデル構築用データの最小値を下回る量の期待値を表します。
  • ここで、 ξ は定数(モデル構築用データの標準偏差の1%ほどが推奨)、Φ 、 ϕ はそれぞれ、標準正規分布の累積分布関数 (CDF) および 確率密度関数 (PDF) を表しています。

以上

3.t分布とEH分布

もくじ

3.1 t分布

1.t分布とは

3.2 EH分布の背景

1.EH分布とは
2.ロバストなベイズ的回帰分析
  • モデルの仮定から大きく逸脱したデータのことを一般的に外れ値と言います．
  • 多くの統計手法は外れ値に大きな影響を受けてしまうことが知られており，推定結果に深刻なバイアスが生じる危険性があります．
  • 簡単な例として以下のような回帰分析の例を考えてみます．(データは擬似的により生成したものです．)
  • 上図では赤線が真の回帰直線を表します．
  • 青線が通常の線形回帰分析によって推定された直線ですが，右上部にある外れ値らしきデータに影響を受けて真の傾きよりも大きい値が推定されてしまっています．
  • この場合の簡単な対処法としては外れ値らしきデータを取り除いてから解析をすることですが，外れ値の特定は必ずしも容易ではありません．
  • 今回のようにデータを図示できる場合は視覚的に外れ値を見つけられる可能性がありますが，説明変数が複数ある一般のケースにおいて外れ値の特定は困難です．
  • そのため，回帰モデルのロバスト(頑健)な推定と外れ値の特定(および除去)を同時に実行するための手法が必要不可欠となります．
3.ベイズ的回帰分析
  • 頻度論的アプローチによるロバストな手法についてはさまざまな研究の蓄積があります．
  • 一方で，ベイズ的な推定を考えた場合，推定結果の不確実性評価が容易であったり，より複雑な統計モデルへの拡張が容易であるなど，さまざまなメリットがあります．そこで，本研究では外れ値がある場合のロバストなベイズ的回帰分析に注目していきます．
  • 以下の(線形)回帰モデルを考えます．

    • $\Large \displaystyle y_i=x_i^t\beta+\sigma \epsilon_i\ \ \epsilon_i\sim N(0,1)$

  • ここで「$y_i$」は応答変数(response variable)，「$x_i$」は共変量(covariate)ベクトルです．
  • また，「$\sigma$」は誤差項の標準偏差です．
  • このモデルの仮定のもと，「$(\beta,\sigma)$」に対する適当な事前分布を導入することで， 「$(\beta,\sigma)$」の事後分布を導出することができます．
  • しかし，誤差項に対する正規分布の仮定は外れ値の影響を受けやすいことが知られています．
  • 広く用いられている改善策としては標準正規分布の代わりにt分布を利用するものがあります．
  • 分布は正規分布よりも裾が重いため外れ値の影響を受けにくくなります．
  • 一方で，t分布は"自由度"と呼ばれるパラメータを持っており，実際に回帰モデルを推定する際にはその自由度パラメータを決定する必要があります．
    • どのような自由度を具体的に設定するべきなのか
    • どのような意味で外れ値の影響を受けにくくなるのか
  • についてはあまり明らかにされておりません．
  • そこで本研究では，ベイズ分析が外れ値に対して頑健であることを定義し，それを達成するために有効な誤差分布をデザインすることを目指しました．

3.3 EH分布は「自由度0のt分布」の誤差分布

1.ベイズ分析におけるロバスト性
  • Dを観測データの集合とし，D*を外れ値を含まないデータの集合と定義します．
  • (もちろんD*がどのような集合になるかについて我々は知ることができません．)
  • また，それぞれのデータに基づく事後分布をそれぞれ「$\pi(\beta,\sigma|D)$」「$\pi(\beta,\sigma|D^*)$」と定義しておきます．
  • 本論文では，ベイズ分析におけるロバスト性として以下のような定義を採用しました．
  • 【外れ値が回帰直線から離れていくにつれ，「$\pi(\beta,\sigma|D)$」が「$\pi(\beta,\sigma|D^*)$」に収束する】
  • これは，どのデータが外れ値かを自動的に検出し，その情報を事後分布から取り除くことが可能であることを示唆しています．
  • すなわち，外れ値でないデータだけを自動的に集めて事後分布を計算できていることを表しています．
2.新しい誤差分布
  • 本研究では新しい誤差分布として以下のような分布を提案しました.

    • $\Large \displaystyle f_{EH}(\epsilon_i)=(1-s)\phi(\epsilon_i;0,1)+s\int_{0}^{\infty}\phi(\epsilon_i;0,u_i)H(u_i,\gamma)du_i$

    $\Large \displaystyle =\frac{\gamma}{1+u_i}\frac{1}{\{1+\log(1+u_i)\}^{1+\gamma}},\ \ \ u_i>0$

  • ただし，「$\phi(\epsilon_i;0,u)$」は、平均「$0$」，分散「$u$」の正規分布の密度関数を表します．
  • また「$s$」は未知の混合パラメータで，「$\gamma$」は形状パラメータです．
  • 上記の分布では，第一項は通常の正規分布になっている一方で，第二項が非常に裾の厚い分布になっており，外れ値の影響を吸収するための分布になっています．
  • 提案分布の具体的な裾の厚さについては，十分大きな に対して

以上

4.深層t過程

もくじ

4.1 t分布とは

1.t分布とは
  • 正規分布は、母平均μと母分散σ2のみで示されたわけですが、一般にこれらの値は知られていない場合が多いようです。
  • そこで、母平均の代わりに標本平均で代用します。
  • ここでは、この標本平均がどれほど母平均に近いかが問題となるわけです。
  • そこで、この標本平均x~はどのような分布を示すのでしょうか？
  • 「$N(\mu,\sigma^2)$」の正規分布を示す母集団から、n 個の標本を取り出しその標本平均を「$\tilde{x}$」とおくと、 「$\tilde{x}$」は「$N(\mu,\frac{\sigma^2}{n})$」の正規分布を示す。
  • つまり、母集団の分布が下図のような正規分布「$N(\mu,\sigma^2$」だとすると、
  • 標本平均の分布は下図のような正規分布「$N(\mu,\frac{\sigma^2}{n})$」となります。
  • つまり、「$\frac{(\tilde{x}-\mu)}{\sqrt{\frac{\sigma^2}{n}}}=(\tilde{x}-\mu)\frac{\sqrt{n}}{\sigma}$」は正規分布「$N(0,1)$」を示します。
  • このように、標本平均が正規分布を示すには標本数が非常に多いということが前提となります。

4.2 t過程回帰モデル

1.t過程回帰モデルとは
  • ガウス過程と同様に、確率過程「$f_{TP}$」が、入力の任意の有限集合「$X$」に対して、 関数値のベクトル「$f_{TP}$」の分布が多変量「$t$」分布であるとき、t過程という。
  • 「$f_{TP}$」の分布は、「$t$」分布であるので、自由度「$\nu\in R^{+}\\[0, 2] $」 ν∈R+[0, 2] と平均関数 µ および、共分散 (カーネル) 関数を要素とする共分散行列 K を指定するによって完全に特定できる。
  • これを、「$f_{TP} \sim T(\nu,\mu(x),K(x,x_0))$」と書く。

以上

5.グラフ描画

「Scilab使用例(4) -グラフ描画の基本-」から引用。

もくじ

5.1 グラフ描画の基礎

1.ウィンドウを1行2列に分割してグラフ描画
  • グラフを作成する為には描画するデータが必要ですが、要素は行方向, 列方向どちらに並べておいても(描画する際に転置すれば良いので)問題ありません。 
    • //ウィンドウを1行2列に分割した左側部分を選択
subplot(1,2,1);

x=[0:0.1:2*%pi]';
plot2d(sin(x));

//右側部分を選択
subplot(1,2,2);

x=[0:0.1:2*%pi]';
plot2d(x,sin(x));
2.ウィンドウを1行4列に分割してグラフ描画
  • plot2d：指定した座標を直線で繋いでいく
  • plot2d2：次の点になるまで同じy座標となる階段状のカーブを描く
  • plot2d3：各点をx軸からの縦線で表現
  • plot2d4：plot2dと同じだが、点を矢印で繋いでいく
  • という風に異なるコマンドを用いることにより見せ方の異なるグラフの描画も可能です。
  • その他にも棒グラフ, 円グラフ, 散布図や等高線図など様々な種類のグラフが用意されていますので、順番に掘り下げた記事をあげていきたいと思います。 
    • clf()
x=[0:0.1:2*%pi]';
subplot(1,4,1); plot2d(x,sin(x));
subplot(1,4,2); plot2d2(x,sin(x));
subplot(1,4,3); plot2d3(x,sin(x));
subplot(1,4,4); plot2d4(x,sin(x));

以上