📖

2023/07/08に公開

# はじめに

## ライブラリのインポート

``````import os
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
import arviz as az

import pymc as pm
import pytensor.tensor as pt

# Initialize random number generator
RANDOM_SEED = 8927
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-darkgrid")
``````

## データの準備

アヒル本８章のデータを使用します。

``````df = pd.read_csv("../input/RStanBook/chap08/input/data-salary-2.txt")
df["KID"] = df["KID"] -1

df_X = df["X"]
df_y = df["Y"]
KID = df["KID"].values
``````

## モデルの定義

### 階層なし

``````with pm.Model() as model_normal:
# coords

# data
X = pm.Data("X", df_X.values, dims=('data'), mutable=True)
y = pm.Data("y", df_y, dims=('data',), mutable=True)

# prior
intercept = pm.Normal("intercept", mu=0, sigma=100)
coef = pm.Normal("coef", mu=0, sigma=100)
sigma_Y = pm.HalfNormal("sigma_Y", sigma=100)

# calc mu
mu = pm.Deterministic("mu", intercept + pt.dot(X, coef), dims=('data',))

# likelihood
obs = pm.Normal("obs", mu=mu, sigma=sigma_Y, observed=y, dims=('data',))
``````
``````pm.model_to_graphviz(model_normal)
``````

### 階層モデル

KIDのグループごとに係数と切片を推論します。ここで、KIDの数だけ`dims=('n_k',)`で次元を指定しています。

``````with pm.Model() as model_hierarchical:
# coords

# data
X = pm.Data("X", df_X.values, dims=('data',), mutable=True)
y = pm.Data("y", df_y, dims=('data',), mutable=True)
kid = pm.Data("kid", KID, dims=('data',), mutable=True)

# prior
intercept_mu_common = pm.Normal("intercept_mu_common", mu=0, sigma=100)
intercept_sigma_common = pm.HalfNormal("intercept_sigma_common", sigma=100)
coef_mu_common = pm.Normal("coef_mu_common", mu=0, sigma=100)
coef_sigma_common = pm.HalfNormal("coef_sigma_common", sigma=100)
sigma_Y = pm.HalfNormal("sigma_Y", sigma=100)

# kidごと
coef = pm.Normal("coef", mu=coef_mu_common, sigma=coef_sigma_common, dims=('n_k',))
intercept = pm.Normal("intercept", mu=intercept_mu_common, sigma=intercept_sigma_common, dims=('n_k',))

# calc mu
mu = pm.Deterministic("mu", intercept[kid] + coef[kid]*X, dims=('data',))

# likelihood
obs = pm.Normal("obs", mu=mu, sigma=sigma_Y, observed=y, dims=('data',))
``````
``````pm.model_to_graphviz(model_hierarchical)
``````

## MCMC

``````with model_normal:
idata_normal = pm.sample(3000,
tune=3000,
nuts_sampler="numpyro",
chains=4,
idata_kwargs={"log_likelihood": True},
random_seed=42,
target_accept=0.90,
return_inferencedata=True)

az.summary(idata_normal, round_to=2)
``````

``````with model_hierarchical:
idata_hierarchical = pm.sample(3000,
tune=3000,
nuts_sampler="numpyro",
chains=4,
idata_kwargs={"log_likelihood": True},
random_seed=42,
target_accept=0.99,
return_inferencedata=True)

az.summary(idata_hierarchical, round_to=2)
``````

## 可視化

``````az.plot_trace(idata_normal);
``````

``````az.plot_trace(idata_hierarchical);
``````

## モデルの比較

### LOO

• rank : モデルのランキング（０が一番いいモデル）
• elpd_loo：looの値
• p_loo：ペナルティ項（雑には有効パラメータ数（汎関数分散）のような値）の値
• elpd_diff：一番いいモデルとの差
• weight : 各モデルに割り当てられた重み, 各モデルが正しい確率として解釈可能
• se : looの標準偏差
• dse : ？？
• warning : 値の信頼性（今回は階層モデルで一部発散しているのでTrueになっている）
• scale : 表示されている値のスケール
``````df_comp_loo = az.compare({"hierarchical": idata_hierarchical, "normal": idata_normal})
df_comp_loo
``````

``````az.plot_compare(df_comp_loo, insample_dev=False);
``````

### WAIC

``````df_comp_waic = az.compare({"hierarchical": idata_hierarchical, "normal": idata_normal}, ic="waic")
df_comp_waic
``````

``````az.plot_compare(df_comp_waic, insample_dev=False);
``````