Python ですらすらわかるベイズ推論「超」入門　メモ

目的

• 統計知識を Python でのコーディングで活かすための実装理解

概要

• ベイズ推論を Python の PyMC というライブラリを駆使して数式をなるべく使わずに実践的に理解していく。

まとめ

第1章 確率分布を理解する

1.1 ベイズ推論における確率分布の必要性

1.2 確率変数と確率分布

1.3 離散分布と連続分布

1.4 PyMCによる確率モデル定義とサンプリング

1.5 サンプリング結果分析

1.6 確率分布とPyMCプログラミングの関係

第2章 よく利用される確率分布

2.1 ベルヌーイ分布(pm.Bernoulliクラス)

• 2 つの事象を 1 回起こした際に得られる分布
  • ex) コインを 1 回投げた場合の裏表の確率

2.2 二項分布(pm.Binomial クラス)

• ベルヌーイ分布を複数回行った場合の拡張
  • ex) コインを複数回投げた場合に表が出る回数の分布

2.3 正規分布(pm.Normal クラス)

• 左右対称の連続分布であるデータの標本数を増加させると基本的に正規分布に近似できる

2.4 一様分布(pm.Uniform クラス)

• 区間内での確率が一定な分布

2.5 ベータ分布(pm.Beta クラス)

• 2項分布の母数が従う分布

2.6 半正規分布(pm.HalfNormal クラス)

• 正規分布を正の値に限定した場合にできる分布
  • ex) 標準偏差の正規分布を求める際には正の値である必要がある

コラム HDI と CI の違い

• HDI: 分布中で N % の面積となる部分を y 軸から囲い込んで x の範囲とする
• CI: 定義域内で N % の面積となる部分を x の定義域から計算し、上下の部分を除外する

第3章 ベイズ推論とは

3.1 ベイズ推論利用の目的

• 事前分布としてある分布を仮定し、実際に起こった事象と重ね合わせることで事後分布として予測分布を導出

3.2 問題設定

3.3 最尤推定による解法

• 最尤推定は実際に起こった事象に対して、最も起こりやすい母数を導出

3.4 ベイズ推論による解法

• 最尤推定とは異なり、母数ではなく分布を導出

3.5 ベイズ推論の精度を上げる方法

• 観測値を増やす
• 事前分布を工夫する

3.6 ベイズ推論の活用例

• ABテスト
• ベイズ回帰モデル
• IRT

第4章 はじめてのベイズ推論実習

4.1 問題設定 (再掲)

• くじ引きで当たる確率を 2 通りの方法で導出

4.2 最尤推定

• Pytorch を用いて対数尤度関数を導出して最大の確率を求める
• 対数尤度関数にマイナスをかけた値を損失として繰り返し処理での微分計算対象としている

import torch # ライブラリインポート

def log_lh(p): # 対数尤度関数
    return (2 * torch.log(p) + 3 * torch.log(1-p))

num_epochs = 40 # 繰り返し回数
lr = 0.01           # 学習率

# パラメータ初期値 (p=0.1)
p = torch.tensor(0.1, dtype=torch.float32, requires_grad=True)

logs = np.zeros((0,3))
for epoch in range(num_epochs):
    loss = -log_lh(p)       # 損失計算
    loss.backward()         # 勾配計算
    with torch.no_grad():
        p -= lr * p.grad    # パラメータ修正
        p.grad.zero_()      # 勾配値の初期化
    log = np.array([epoch, p.item(), loss.item()]).reshape(1,-1)
    logs = np.vstack([logs, log])

plt.rcParams['figure.figsize'] = (8, 4)
fig, axes = plt.subplots(1, 2)
axes[0].plot(logs[:,0], logs[:,1])
axes[0].set_title('p(確率値)')
axes[1].plot(logs[:,0], logs[:,2])
axes[1].set_title('loss(損失)')
plt.tight_layout();

4.3 ベイズ推論 (確率モデル定義)

# コンテキスト定義
model1 = pm.Model()

with model1:
    # pm.Uniform: 一様分布
    p = pm.Uniform('p', lower=0.0, upper=1.0)

    # pm.Bernoulli: ベルヌーイ分布
    X_obs = pm.Bernoulli('X_obs', p=p, observed=X)

g = pm.model_to_graphviz(model1) # 可視化
display(g)

4.4 ベイズ推論 (サンプリング)

• モデル定義したコンテキスト内でサンプリング
• 定常状態の確率変数だけを取り出したいので、最初の方のサンプルは除外

with model1:
    idata1_1 = pm.sample(
        # 乱数系列の数(デフォルト2)
        chains=3,
        # 捨てるサンプル数(デフォルト1000)
        tune=2000,
        # 取得するサンプル数(デフォルト1000)
        draws=2000,
        random_seed=42)

4.5 ベイズ推論 (結果分析)

• サンプリングした各データセットに対して横軸：確率変数、縦軸：頻度としてグラフ描画で結果を確認可能
• 確率の収束過程もグラフ描画で追うことが可能

axes = az.plot_trace(idata1_1, compact=False)
plt.tight_layout();

axes = az.plot_trace(idata1_2, compact=False)
plt.tight_layout(); # idata_1_2 は 2 群のデータセットなので、グラフは 2 本

• 事後分布を plot_posterior 関数で描画可能

plt.rcParams['figure.figsize'] = (6, 6)
ax = az.plot_posterior(idata1_2)
ax.set_xlim(0, 1)
ax.set_title('ベイズ推論結果　初期版');

• summary 関数で各統計値を表形式で表示可能

mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
p	0.377	0.162	0.083	0.678	0.006	0.004	844.000	960.000

4.6 ベイズ推論 (二項分布バージョン)

• ベルヌーイ分布と同様に PyMC の関数を用いてベイズ推論
• モデル定義 (可視化) → サンプリング → plot_trace でサンプル数による確率の収束観察 → 事後分布描画 → summary 関数で各統計値のまとめ
• ベルヌーイ分布とは異なり、各状態が観測された回数のみを記述すれば良いのでモデル定義のところだけ異なる

# コンテキスト定義
model2 = pm.Model()

with model2:
    # pm.Uniform: 一様分布
    p = pm.Uniform('p', lower=0.0, upper=1.0)

    # pm.Binomial: 二項分布
    # p: 成功確率
    # n: 試行数
    X_obs = pm.Binomial('X_obs', p=p, n=5, observed=2) # ここで観測された事象をそのまま書ける

4.7 ベイズ推論 (試行回数を増やす)

• 試行回数を増やすと確率が同様に確からしい場合の確率に近い値に収束する可能性が高まる

# コンテキスト定義
model3 = pm.Model()

with model3:
    # pm.Uniform: 一様分布
    p = pm.Uniform('p', lower=0.0, upper=1.0)

    # pm.Binomial:  二項分布
    # p: 成功確率
    # n: 試行数
    X_obs = pm.Binomial('X_obs', p=p, n=50, observed=20)

• 試行回数 5 回

• 試行回数 50 回

↑確かに左のグラフでは試行回数 50 回の方が観測された事象から予測される確率 0.4 周辺にピークが存在

右の

グラフでも試行回数 50 回の方が 0.4 に近い値に収束している (ばらつきが小さい)

• summary 関数の結果からも HDI の範囲を狭められている

mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
p	0.403	0.068	0.277	0.535	0.002	0.002	854.000	1163.000

4.8 ベイズ推論 (事前分布の変更)

• 事前分布を狭めることでも HDI の範囲を狭めることができる

# コンテキスト定義
model4 = pm.Model()

with model4:
    # 確率モデル定義

    # 一様分布のパラメータを変更
    p = pm.Uniform('p', lower=0.1, upper=0.9)

    # 5回中2回あたりという観測値はそのまま
    X_obs = pm.Binomial('X_obs', p=p, n=5, observed=2)

    # サンプル値取得
    idata4 = pm.sample(random_seed=42)

mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
p	0.442	0.169	0.134	0.739	0.006	0.004	783.000	976.000

4.9 ベータ分布で直接確率分布を求める

• 事前分布：一様分布、観測事象：二項分布の事後分布はベータ分布であるが、これをベイズ推論を使わずにそのまま求めてみる

# 真のベータ関数の定義
from scipy import stats
alpha = 20 + 1
beta = 30 + 1
true_beta = stats.beta(alpha, beta)

plt.rcParams['figure.figsize'] = (6, 6)
# ベイズ推論結果の可視化
# idata3は4.7節で計算した結果を利用
ax = az.plot_posterior(idata3)
ax.lines[0].set_label('ベイズ推論結果')

# 真のベータ関数の可視化
x = np.linspace(*ax.get_xlim())
ax.plot(x, true_beta.pdf(x), color='orange', label='真値')
ax.legend(loc='center right');

第5章 ベイズ推論プログラミング

• 実際のデータセットでどのようにベイズ推論を適用できるかを検証

5.1 データ分布のベイズ推論

• seaborn 経由で iris データセットをインポート

# アイリスデータセットの読み込み
df = sns.load_dataset('iris')

# 先頭5件の確認
display(df.head())

#  speciesの分布確認
df['species'].value_counts()

	sepal_length	sepal_width	petal_length	petal_width	species
0	5.100	3.500	1.400	0.200	setosa
1	4.900	3.000	1.400	0.200	setosa
2	4.700	3.200	1.300	0.200	setosa
3	4.600	3.100	1.500	0.200	setosa
4	5.000	3.600	1.400	0.200	setosa

• ヒストグラムの描画
  • ヒストグラムのカーネル密度推定 (KDE) の関数とベイズ推論の分布がどの程度近くなるかどうかが今回のゴール

# setosaの行のみ抽出
df1 = df.query('species == "setosa"')

bins = np.arange(4.0, 6.2, 0.2)
# ヒストグラムを描画
sns.histplot(df1, x='sepal_length', bins=bins, kde=True)
plt.xticks(bins);

• 変数 X と値自体の確認

# sepal_length列の抽出
s1 = df1['sepal_length']

# NumPy変数の1次元配列に変換
X = s1.values

# 統計情報の確認
print(s1.describe())

# 値の確認
print(X)

count   50.000
mean     5.006
std      0.352
min      4.300
25%      4.800
50%      5.000
75%      5.200
max      5.800
Name: sepal_length, dtype: float64
[5.100 4.900 4.700 4.600 5.000 5.400 4.600 5.000 4.400 4.900 5.400 4.800
 4.800 4.300 5.800 5.700 5.400 5.100 5.700 5.100 5.400 5.100 4.600 5.100
 4.800 5.000 5.000 5.200 5.200 4.700 4.800 5.400 5.200 5.500 4.900 5.000
 5.500 4.900 4.400 5.100 5.000 4.500 4.400 5.000 5.100 4.800 5.100 4.600
 5.300 5.000]

• sepal_length が正規分布に従うと仮定しているので、正規分布のパラメータ $\mu, \sigma$ を設定しこれを sepal_length の値から推定していく
  • 正規分布の事前分布が設定されていないが、正規分布の事前分布は正規分布であると仮定すれば事後分布も正規分布となるということがわかっているので、パラメータをそれぞれ正規分布に従うものとして事前分布を設定
  • $\sigma$ は正の値である必要があるので、半正規分布に従うこととしている

model1 = pm.Model() # 確率モデルの定義

with model1:
  mu = pm.Normal('mu',mu=0.0, sigma=10.0)
  sigma = pm.HalfNormal('sigma', sigma=10.0)
  X_obs = pm.Normal('X_obs', mu=mu, sigma=sigma, observed=X)

• パラメータのサンプリング
  • サンプリング回数によって値にムラがないのでうまくサンプリングできていることがわかる

• 事後分布の描画

• summary 関数で分析結果の表示
  • mcse_mean, ess_bulk の値が小さいことからもうまくベイズ推論できていることがわかる

mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
mu	5.005	0.050	4.908	5.093	0.001	0.001	1744.000	1399.000
sigma	0.362	0.038	0.294	0.438	0.001	0.001	1519.000	1011.000

• 推論したパラメータを用いてグラフ描画

def norm(x, mu, sigma): # 正規分布関数
    y = (x-mu)/sigma
    a = np.exp(-(y**2)/2)
    b = np.sqrt(2*np.pi)*sigma
    return a/b

x_min = X.min()
x_max = X.max()
x_list = np.arange(x_min, x_max, 0.01)
y_list = norm(x_list, mu_mean1, sigma_mean1) # 推論結果から得られる分布の平均値をパラメータとして用いて分布を描画

delta = 0.2
bins=np.arange(4.0, 6.0, delta)
fig, ax = plt.subplots()
sns.histplot(df1, ax=ax, x='sepal_length',
    bins=bins, kde=True, stat='probability')
ax.get_lines()[0].set_label('KDE曲線')
ax.set_xticks(bins)
ax.plot(x_list, y_list*delta, c='b', label='ベイズ推論結果')
ax.set_title('ベイズ推論結果とKDE曲線の比較')
plt.legend();

• 少ないデータ数でのベイズ推論

  • データ数が少なくともある程度の精度でベイズ推論可能
  • データ数を 5 個とした場合の各パラメータの分布

  

  • 平均値はほとんど変わらないものの、HDI の範囲がかなり増えている (データ数を減らした影響)

• $\sigma$ を $\tau = 1 / \sigma^2$ として正規分布のサンプリングを行うことも可能 (正規分布の式中では $1/\sigma^2$ の形で使われるため)

5.2 線形回帰のベイズ推論

• アイリスデータセットで花の種類を特定した場合、２つの項目値には正の相関があり。線形回帰(単回帰)に従うと見なせる。 当実習では、上の点を前提とした上で、最適な回帰式をベイズ推論で求める。
• $Y = \alpha X_n + \beta + \epsilon_n$ として線形回帰のパラメータをベイズ推論で決めていく

model1 = pm.Model()

with model1:
    # 確率変数alpha、betaの定義(一次関数の傾きと切片)
    alpha = pm.Normal('alpha', mu=0.0, sigma=10.0)
    beta = pm.Normal('beta', mu=0.0, sigma=10.0)

    # 平均値muの計算
    mu = alpha * X + beta

    #　誤差を示す確率変数epsilonの定義
    epsilon = pm.HalfNormal('epsilon', sigma=1.0)

    # 観測値を持つ確率変数はY_obsとして定義
    Y_obs= pm.Normal('Y_obs', mu=mu, sigma=epsilon, observed=Y)

• ↑ だと $X$ が図中に出てこないのが問題なので、定数として $X, Y$ を定義

model2 = pm.Model()

with model2:
    # X, Yの観測値をConstantDataとして定義
    X_data = pm.ConstantData('X_data', X)
    Y_data = pm.ConstantData('Y_data', Y)

    # 確率変数alpha、betaの定義(一次関数の傾きと切片)
    alpha = pm.Normal('alpha', mu=0.0, sigma=10.0)
    beta = pm.Normal('beta', mu=0.0, sigma=10.0)

    # 平均値muの計算
    mu = pm.Deterministic('mu', alpha * X_data + beta)

    #　誤差を示す確率変数epsilonの定義
    epsilon = pm.HalfNormal('epsilon', sigma=1.0)

    # 観測値を持つ確率変数はobsとして定義
    obs = pm.Normal('obs', mu=mu, sigma=epsilon, observed=Y_data)

• それぞれのパラメータの分布分析

• 事後分布描画

• 推論結果から得られるパラメータを用いて線形回帰のグラフ描画

for y_pred in y_preds:
    plt.plot(x_values, y_pred, lw=1, alpha=0.01, c='g')
plt.scatter(X, Y)
plt.title('ベイズ推論による回帰直線')
plt.xlabel('sepal_length')
plt.ylabel('sepal_width');

• 少ないデータ数でもある程度の精度のベイズ推論が可能
  • 2 つのサンプルの概形がほとんど同じである一方で、誤差の値が大きな値になってしまっている

• target_accept の値をデフォルトの値 (0.80) にすると、推論中に発散してしまいうまく推論できない
• 一方で target_accept の値を 1 に近い値にすると条件が厳しいので精度が上がる一方で、サンプリングに時間がかかる

5.3 階層ベイズモデル

• アイリスデータセットで３種類の花の回帰直線を同時に求める。データ数が少ないという条件の下で、３つの回帰直線は共通の傾向がある前提とする。この場合、階層ベイズモデルの問題に帰着する。
• 1 層のベイズモデルではサンプル数が少ない場合に精度よく推論するのが難しい → ターゲットとなるデータセット以外に類似するデータセットが存在する場合に推論のデータセットとして含めることで精度の良い推論を実現 →  事前分布のパラメータをさらに事前分布として推論をかけ、多層構造とすることで推論精度を上げる (階層ベイズモデル)
• 1 種類の線形回帰の推論の問題は、傾きと切片の値をある固定の正規分布を事前分布として推論していた
• 今回はそれは異なり、種類ごとに事前分布のパラメータ自体を変更する必要があるので、事前分布のパラメータ自体にも正規分布からサンプリングを行ってパラメータ決定している

model1 = pm.Model()

with model1:
    # X, Yの観測値をConstantDataとして定義(通常ベイズと共通)
    X_data = pm.ConstantData('X_data', X)
    Y_data = pm.ConstantData('Y_data', Y)

    # クラス変数定義(階層ベイズ固有)
    cl_data = pm.ConstantData('cl_data', cl)

    # 確率変数alphaの定義(階層ベイズ固有)
    a_mu = pm.Normal('a_mu', mu=0.0, sigma=10.0)
    a_sigma = pm.HalfNormal('a_sigma',sigma=10.0)
    alpha = pm.Normal('alpha', mu=a_mu, sigma=a_sigma, shape=(3,))

    # 確率変数betaの定義(階層ベイズ固有)
    b_mu = pm.Normal('b_mu', mu=0.0, sigma=10.0)
    b_sigma = pm.HalfNormal('b_sigma', sigma=10.0)
    beta = pm.Normal('beta', mu=b_mu, sigma=b_sigma, shape=(3,))

    #　誤差epsilon(通常ベイスと共通)
    epsilon = pm.HalfNormal('epsilon', sigma=1.0)

    # muの値は、cl_dataによりindexを切り替えて計算(階層ベイズ固有)
    mu = pm.Deterministic('mu', X_data * alpha[cl_data] + beta[cl_data])

    # mu, epsilonを使って観測値用の確率モデルを定義(通常ベイスと共通)
    obs = pm.Normal('obs', mu=mu, sigma=epsilon, observed=Y_data)

• 推論結果確認

            mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
  

alpha[0] 0.619 0.748 -0.488 2.061 0.041 0.029 451.000 349.000 1.010
alpha[1] 0.435 0.228 0.042 0.878 0.013 0.009 305.000 474.000 1.020
alpha[2] 0.271 0.205 -0.138 0.650 0.009 0.007 540.000 420.000 1.000
beta[0] 0.212 3.714 -6.136 6.360 0.202 0.143 456.000 351.000 1.010
beta[1] 0.307 1.292 -2.155 2.629 0.074 0.052 301.000 486.000 1.020
beta[2] 1.110 1.357 -1.515 3.687 0.059 0.042 541.000 418.000 1.000

• 抽出サンプルが耀っているにも関わらず、青線は他の線とほぼ平行になっている

• PyMC でどの要素まで表示するか？ → 慣れるまではなるべくすべての計算過程を表示していた方が良い (右より左)

  

  

5.4 潜在変数モデル

• 現実のデータでは、異なる種別のデータが混合していることが多く、それらのデータセットを潜在変数モデルを用いることで分けることが可能

# 変数の初期設定

# 何種類の正規分布モデルがあるか
n_components = 2

# 観測データ件数
N = X.shape

model1 = pm.Model()

with model1:
    # Xの観測値をConstantDataとして定義
    X_data = pm.ConstantData('X_data', X)

    # p: 潜在変数が1の値をとる確率
    p = pm.Uniform('p', lower=0.0, upper=1.0)

    # s: 潜在変数pの確率値をもとに0, 1のいずれかの値を返す
    s = pm.Bernoulli('s', p=p, shape=N)

    # mus: ２つの花の種類毎の平均値
    mus = pm.Normal('mus', mu=0.0, sigma=10.0, shape=n_components)

    # taus: 2つの花の種類毎のバラツキ
    # 標準偏差sigmasとの間にはtaus = 1/(sigmas*sigmas)の関係がある
    taus = pm.HalfNormal('taus', sigma=10.0, shape=n_components)

    # グラフ描画など分析でsigmasが必要なため、tausからsigmasを求めておく
    sigmas = pm.Deterministic('sigmas', 1/pm.math.sqrt(taus))

    # 各観測値ごとに潜在変数からmuとtauを求める
    mu = pm.Deterministic('mu', mus[s])
    tau = pm.Deterministic('tau', taus[s])

    # 正規分布に従う確率変数X_obsの定義
    X_obs = pm.Normal('X_obs', mu=mu, tau=tau, observed=X_data)

• これまでの推論モデルは固定のインデックス変数を用いていたが、潜在変数モデルでは可変の潜在変数の値によってインデックス変数を可変にしている
• サンプリング結果
  • 潜在変数 p：0.5 あたりでピークになっているので妥当そう
  • $\mu , \sigma$ ：2 つのピークの部分がきれいに分かれていて妥当そう

• 事後分布の描画

• 各統計値

  mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
  p 0.483 0.110 0.254 0.662 0.011 0.008 100.000 178.000 NaN
  mus[0] 1.357 0.060 1.255 1.482 0.005 0.003 176.000 330.000 NaN
  mus[1] 2.023 0.095 1.848 2.203 0.010 0.007 101.000 131.000 NaN
  sigmas[0] 0.242 0.036 0.181 0.307 0.002 0.002 274.000 385.000 NaN
  sigmas[1] 0.298 0.055 0.205 0.400 0.005 0.004 109.000 177.000 NaN

• 元のヒストグラムの推測関数とのグラフの比較

  • KDE は花の種類がわかった上で近似曲線を描いているので、それとほとんど一致しているということはベイズ推論が正しくできているということ

  

• ↑のモデルで推論に使う変数を $\sigma$ にするとどうなるのか？

  

• 潜在変数 p がサンプリング’の初期でほとんど 1 になっており、うまくグループ分けできていない

• 変数の順序関係を固定するとラベルスイッチが起きない

  # 変数の初期設定
  
  # 何種類の正規分布モデルがあるか
  n_components = 2
  
  # 観測データ件数
  N = X.shape
  
  model3 = pm.Model()
  
  with model3:
      # Xの観測値をConstantDataとして定義
      X_data = pm.ConstantData('X_data', X)
  
      # p: 潜在変数が1の値をとる確率
      p = pm.Uniform('p', lower=0.0, upper=1.0)
  
      # s: 潜在変数　pの確率値をもとに0, 1のいずれかの値を返す
      s = pm.Bernoulli('s', p=p, shape=N)
  
      # mus: ２つの花の種類毎の平均値
      mu0 = pm.HalfNormal('mu0', sigma=10.0)
      delta0 = pm.HalfNormal('delta0', sigma=10.0)
      mu1 = pm.Deterministic('mu1', mu0+delta0)
      mus = pm.Deterministic('mus',pm.math.stack([mu0, mu1]))
  
      # taus: 2つの花の種類毎のバラツキ
      # 標準偏差sigmasとの間にはtaus = 1/(sigmas*sigmas)の関係がある
      taus = pm.HalfNormal('taus', sigma=10.0, shape=n_components)
  
      # グラフ描画など分析でsigmasが必要なため、tausからsigmasを求めておく
      sigmas = pm.Deterministic('sigmas', 1/pm.math.sqrt(taus))
  
      # 各観測値ごとに潜在変数からmuとtauを求める
      mu = pm.Deterministic('mu', mus[s])
      tau = pm.Deterministic('tau', taus[s])
  
      # 正規分布に従う確率変数X_obsの定義
      X_obs = pm.Normal('X_obs', mu=mu, tau=tau, observed=X_data)
  
  g = pm.model_to_graphviz(model3)
  display(g)
  

第6章 ベイズ推論の業務活用事例

6.1 ABテストの効果検証

• 問題設定
  • Web ページの AB テストの効果検証
• Web ページなので、表示された回数に対するクリック率を比較したい
• 各バージョンのクリック率は2項分布に従うものの、サンプル数が少なすぎるので、ベイズ推定で分布を求める
• それぞれの差分の分布も導出することで一方が他方よりも優れているかどうかを判定

model_s = pm.Model()

with pm.Model() as model_s:
    # 事前分布として一様分布を採用
    p_s_a = pm.Uniform('p_s_a', lower=0.0, upper=1.0)
    p_s_b = pm.Uniform('p_s_b', lower=0.0, upper=1.0)

    # 二項分布で確率モデルを定義
    # n:表示数　observed:ヒット数　とすればよい
    obs_s_a = pm.Binomial('obs_s_a', p=p_s_a, n=40, observed=2)
    obs_s_b = pm.Binomial('obs_s_b', p=p_s_b, n=25, observed=2)

    # 新たな確率変数として二つの確率変数の差を定義
    delta_prob_s = pm.Deterministic('delta_prob_s', p_s_b - p_s_a)

# 確率モデル構造可視化
g = pm.model_to_graphviz(model_s)
display(g)

• それぞれの変数が2群のサンプリングであまりぶれていないので妥当そう

• AB の差分のうち、マイナスになる確率が 28.60 % であり、ページを変えたことによる効果は必然とは言い切れない

• 異なる場合では以下のように有意に一方の方が’効果がありそう

  

• 今回の場合だと、事前分布が一様分布で観測対象のモデルが2項分布なので事後分布はベータ分布に従うと考えられる

  # 画面A 成功2回　失敗3８回
  alpha_a = 2 + 1
  beta_a = 38 + 1
  
  # 画面B 成功2回　失敗23回
  alpha_b = 2 + 1
  beta_b = 23 + 1
  
  model_s2 = pm.Model()
  with model_s2:
      # 確率モデル定義
      # pm.Beta: ベータ分布
      # alpha: 注目している試行の成功数+1
      # beta: 注目している試行の失敗数+1
      p_a = pm.Beta('p_a', alpha=alpha_a, beta=beta_a)
      p_b = pm.Beta('p_b', alpha=alpha_b, beta=beta_b)
  
      # サンプル値取得
      samples_s2 = pm.sample_prior_predictive(random_seed=42, samples=10000)
  
  # サンプル値抽出
  p_a_samples_s2 = samples_s2['prior']['p_a'].values.reshape(-1)
  p_b_samples_s2 = samples_s2['prior']['p_b'].values.reshape(-1)
  delta_a_b_s2 = p_b_samples_s2 - p_a_samples_s2
  
  # delta_probの値がマイナスであった件数
  n1_s2 = (delta_a_b_s2 < 0).sum()
  
  # 全体サンプル数
  n_s2 = len(delta_a_b_s2)
  
  # 比率計算
  n1_rate_s2 = n1_s2/n_s2
  
  # 可視化
  ax = az.plot_dist(delta_a_b_s2)
  xx, yy = ax.get_lines()[0].get_data()
  ax.fill_between(xx[xx<0], yy[xx<0])
  
  # グラフタイトル
  title = f'鈴木さんケース　画面Aの方がクリック率が高い確率(別解):\
  {n1_rate_s2*100:.02f}%'
  ax.set_title(title, fontsize=12);
  

6.2 ベイズ回帰モデルによる効果検証

• 難聴が子供の音声言語の発達に影響を与える要因を調査したデータを使用

• 目的変数 score

  

• 目的変数と説明変数にデータセットを分離したり、変数の正規化などのデータ整形を実行

• 確率モデルの定義

  # 説明変数リストをpredictorsとして定義
  model1 = pm.Model(coords={'predictors': columns})
  
  with model1:
      # Xは従来のベクトルが行列になる。転置していることに注意
      X_data = pm.ConstantData('X_data', X.T)
  
      # yが回帰モデルの目的変数
      y_data = pm.ConstantData('y_data', y)
  
      # 単回帰のときスカラーだったalphaは重回帰でベクトルになる
      # 要素数はpredictorsにより間接的に指定できる(上でcoordsパラメータを指定した効果)
      alpha = pm.Normal('alpha', mu=0.0, sigma=10.0, dims='predictors')
  
      # betaとepsilonは単回帰の時と同じ(パラメータ値の選定理由は本文で説明)
      beta = pm.Normal('beta', mu=100.0, sigma=25.0)
      epsilon = pm.HalfNormal('epsilon', sigma=25.0)
  
      # muの計算では、単回帰のときのかけ算が内積に変わっている
      mu = pm.Deterministic('mu', alpha @ X_data + beta)
  
      # 正規分布の定義は5.2節の単回帰と同じ
      obs = pm.Normal('obs', mu=mu, sigma=epsilon, observed=y_data)
  
  g = pm.model_to_graphviz(model1)
  display(g)
  

  

• サンプリング

  • alpha に関しては変数の数だけ存在 (上の図からもわかる通り)

• 各統計情報
  • r_hat は 1 に近いほど収束が良いという意味

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
alpha[male]	0.989	1.791	-2.379	4.317	0.036	0.039	2486.000	1298.000	1.000
alpha[siblings]	-1.972	1.818	-5.324	1.408	0.039	0.031	2165.000	1421.000	1.000
alpha[family_inv]	-9.226	2.186	-13.879	-5.537	0.051	0.036	1805.000	1317.000	1.000
alpha[non_english]	-3.765	1.845	-7.349	-0.460	0.038	0.029	2273.000	1568.000	1.000
alpha[prev_disab]	-5.134	1.913	-8.950	-1.683	0.039	0.029	2325.000	1268.000	1.000
alpha[age_test]	1.642	1.852	-1.812	5.026	0.039	0.031	2208.000	1628.000	1.000
alpha[non_severe_hl]	2.854	1.795	-0.674	6.111	0.040	0.033	2014.000	1298.000	1.000
alpha[mother_hs]	0.298	2.023	-3.353	4.288	0.043	0.043	2222.000	1581.000	1.000
alpha[early_ident]	4.193	1.737	0.815	7.435	0.036	0.026	2365.000	1680.000	1.000
alpha[non_white]	-2.848	1.959	-6.647	0.770	0.043	0.032	2083.000	1602.000	1.000

6.3 IRT (Item Response Theory)によるテスト結果評価

• 以下の表形式のデータから各設問の正答率と難易度・受験者の学力を評価したい

• 複数の説明変数に対して一つの目的変数が関係している状況で説明変数を同時に推定する必要がある場合

  	Q001	Q002	Q003	Q004	Q005	Q006	Q007	Q008	Q009	Q010	Q011	Q012	Q013	Q014	Q015	Q016	Q017	Q018	Q019	Q020	Q021	Q022	Q023	Q024	Q025	Q026	Q027	Q028	Q029	Q030	Q031	Q032	Q033	Q034	Q035	Q036	Q037	Q038	Q039	Q040	Q041	Q042	Q043	Q044	Q045	Q046	Q047	Q048	Q049	Q050
  USER0001	0	1	1	1	0	1	1	0	0	0	1	1	1	0	1	1	0	1	0	1	1	0	1	1	0	0	1	1	1	1	0	1	1	1	1	0	1	1	1	0	0	0	1	1	0	1	0	1	1	1
  USER0002	1	0	1	1	1	0	1	1	0	0	1	1	1	0	0	0	1	1	1	0	0	1	1	1	0	1	1	1	1	1	0	0	0	1	1	0	0	1	1	0	0	1	1	0	1	0	0	1	1	0
  USER0003	1	0	1	1	1	1	1	1	0	0	1	1	1	0	1	1	1	1	1	1	1	1	0	1	0	0	1	1	1	1	1	0	0	0	1	1	1	0	1	1	1	1	1	1	1	1	1	1	1	1
  USER0004	1	1	1	1	1	0	1	0	1	0	0	1	1	1	0	1	0	1	0	1	0	0	0	0	0	0	1	0	1	0	1	0	0	0	0	1	1	0	1	0	0	0	1	1	1	1	0	0	0	1
  USER0005	0	1	0	1	0	0	1	0	1	1	0	0	1	0	0	0	1	1	1	0	1	1	0	0	0	0	1	0	0	1	0	0	0	1	1	1	1	1	1	1	1	1	1	0	1	0	0	1	0	0

• 横持ちのデータを縦持ちに変換＆カテゴリカルデータを数値変換し、データ整形

• 目的変数が正答かどうかなので、ロジスティック回帰

• ロジスティック回帰の際の変数 $x$ が問題難易度や回答者の学力・識別力に依存していることを考慮

# 配列の項目定義(ユーザー軸と問題軸の2軸)
coords = {'user': users, 'question': questions}

# 確率モデルインスタンスの定義
model1 = pm.Model(coords=coords)

with model1:
    # 観測値の配列(1:正答　0:誤答)
    response_data = pm.ConstantData('response_data', response)

    # 能力値（受験者ごと）
    theta = pm.Normal('theta', mu=0.0, sigma=1.0, dims='user')

    # 識別力（設問ごと）
    a = pm.HalfNormal('a', sigma=1.0, dims='question')
    # 困難度（設問ごと）
    b = pm.Normal('b', mu=0.0, sigma=1.0, dims='question')

    # logit_pの計算 (2パラメータ・ロジスティックモデル（2PLM）)
    logit_p = pm.Deterministic(
        'logit_p', a[question_idx] * (theta[user_idx] - b[question_idx]))

    # ベルヌーイ分布の定義(1:正答　0:誤答)
    obs = pm.Bernoulli('obs', logit_p=logit_p, observed=response_data)

g = pm.model_to_graphviz(model1)
display(g)

• 推論結果(能力値)と偏差値との関係

  

• 同一の偏差値となっている受験者間の能力値の違いはどこからきているのか

  • 正答している問題が片方の受験者の方が難易度が高かったものの、正答数がおいなじであったために偏差値が同じであった

    w1のshape: (2, 50)
    正解数:
    USER0463    31
    USER0064    31
    dtype: int64
    [[-1.173 -1.419 -0.700 -1.789 -2.339 -1.405 -2.134 -1.840  1.244  1.280
    -0.565 -1.939 -0.674  0.278 -0.086  0.329  0.071 -2.377 -1.763 -1.030
    -1.568 -1.856  1.347 -1.515  1.289  1.071 -1.677 -0.288 -1.937 -2.184
    -0.175  0.013 -0.904 -0.646 -1.263  0.265 -1.546 -1.401 -2.087  0.600
    -0.069 -1.309 -2.847 -1.367  0.616  0.459  0.529 -1.775 -2.003 -1.863]]
    b_meanのshape: (1, 50)
    
    USER0463   -1.204
    USER0064   -1.113
    dtype: float64
    

    

参考サイト

• 書籍情報

https://www.kspub.co.jp/book/detail/5337639.html

• サンプルコード

https://github.com/makaishi2/python_bayes_intro