【改】第４章「聞きづらい行為を他のみんなはどれくらいしているか」のベイズモデリングをPyMC Ver.5 で

この記事は、以前投稿した『 第４章「聞きづらい行為を他のみんなはどれくらいしているか」のPyMC記事 』を改造する記事です。

以前の記事はこちらです。

以前の記事執筆時点では、事後分布を収束させることができず、苦い思いをしました。
しかし、事態は一変します。
生成AIの進化のおかげで、解決策が見つかったのです！

それではテイクツー、スタートです！

リスタート

この記事は、テキスト「たのしいベイズモデリング２」の第４章「聞きづらい行為を他のみんなはどれくらいしているか」のベイズモデルを用いて、PyMC Ver.5で「実験的」に実装する様子を描いた統計ドキュメンタリーです。

テキストは第２章から第４章まで「間接質問法」の各種手法で「違法行為等の経験」の調査データを用いたベイズモデリングに取り組んでいます。
この章では、「Aggregated Response法」(AR法)を用いてデリケートな行為の経験を推論します。

AR法には前書「たのしいベイズモデル」第１章で出会いました。
この当時は自己流モデルの出力がテキストの結果と相違して、苦い思いをしました。
さあ、心新たにPyMCのベイズモデリングを楽しみましょう。

テキストの紹介、引用表記、シリーズまえがき、PyMC等のバージョン情報は、このリンクの記事をご参照ください。

テキストで使用するデータは、R・Stan等のサンプルスクリプトとともに、出版社サイトからダウンロードして取得できます。

サマリー

テキストの概要

執筆者　　　：　秋山隆 先生、豊田秀樹 先生、久保沙織 先生、岡律子 先生
モデル難易度：　★★★・・　（ふつう）

自己評価

評点

$$
\begin{array}{c:c:c}
実装精度 & ★★★★★ & GoooD! \\
結果再現度 & ★★★★★ & 最高！ \\
楽しさ & ★★★★★ & 楽しい！ \\
\end{array}
$$

評価ポイント

• 以前の記事ではポアソン分布モデルのカスタム分布関数で発生するエラーを解消できず、不完全燃焼のまま終演を迎えてしまいました。

• しかし、エラーメッセージとモデルのコードを「ポン！」と生成AIに投げたところ、いとも簡単に動くコードが生成されたのです。
  頼もしい相手に気軽に相談できる時代になりました。

工夫・喜び・反省

• 積み残した宿題を無事にクローズできました！
  やったね✨️

モデルの概要

テキストの調査・実験の概要

■ 間接質問法
間接質問法は、答えにくい質問の回答を得る場合に、回答者のプライバシーを守りつつ、分析に必要なデータを取得できる質問方法です。
テキスト第２章の定義をお借りします。

この章の調査では、大学生のデリケートな行為の経験回数・人数（カウントデータ）を取り扱っています。
上記定義に当てはめると、「調査者が回答を見ても大学生個人の回答値が分からないように配慮しつつ、ある大学生グループ全体の推定経験量を得る方法」と言い換えられます。

■ AR法
間接質問法の手法の１つです。
下図をご覧ください。
質問に対して「真の値 ± マスクの値」で回答します。

携帯番号下１桁のように、真の値に加減算する数値を「マスク」と呼びます。
マスクは「秘匿（マスク）」となり、プライバシー保護として機能するので、回答者は回答しやすくなります。

■ 調査の概要
「性行為の経験人数」（以下、Sと略します）と「カンニングの経験回数」（以下、Cと略します）を調査項目としています。
Google フォームで質問リストを作成して、調査参加者からWeb上で回答を収集し、Sは130件、Cは125件のデータを用いて分析を進めたそうです。

テキストのモデリング

この章ではモデルの掲載はなく、「たのしいベイズモデリング」第１章を参照するように書かれています。
という事で、「たのしいベイズモデリング」第１章のブログもご参照ください！
「テキストのモデリング」の項に、マスクがモデリングにどのように関わるかを詳しく書いています。

以下のモデリング内容は、「たのしいベイズモデリング」第１章を参考にして書きます。

■目的変数と関心のあるパラメータ
目的変数はマスクされた回答値$${y}$$です。
関心のあるパラメータはマスクを除いた真の人数・回数$${\lambda}$$です。

■ テキストのベイズモデル
テキスト「たのしいベイズモデリング」第１章のモデルを引用いたします。
なお、この章では回答値が平均パラメータ$${\lambda}$$のポアソン分布に従うと仮定されている点を、次の数式に反映しています。

■マスク$${m_k}$$、集合$${A}$$、関数$${N(k \subset A)}$$
マスク$${\boldsymbol{m_k}}$$の取りうる値の集合$${\boldsymbol{A}}$$は、回答$${\boldsymbol{y}}$$の値によって変動します。

たとえば回答値$${y=2}$$の場合、加減算する携帯番号下１桁（＝マスク）は、コイン裏＆下１桁９の$${-9}$$から、コイン表＆下１桁２の$${+2}$$までが取りうる値になります。
ではコイン表＆下１桁３の場合、どうなるでしょう。
経験人数・回数の最小値が$${0}$$のため回答は$${3}$$以上になってしまい、回答値$${y=2}$$に収まらないのです。
したがって、マスクの値は$${+3}$$以上を取り得ないのです。
まとめますと、回答値$${y=2}$$の場合、マスクの値の集合$${A=\{ -9, -8, \cdots, -1, -0, +0, 1, 2\}}$$となり、要素数$${N(k \subset A)}$$は$${13}$$となります。

■混合分布であること
$${f(y-m_k \mid \cdot \ )}$$は、$${y}$$の値によって変動する集合$${A}$$の要素の数（＝マスクの数）だけ、確率分布が存在することを示しています。
回答$${\boldsymbol{y}}$$の値によって、混合する「確率分布の数が変動する」のです！
たとえば$${y=2}$$の場合、混合分布は、以下の$${13}$$個のポアソン分布の確率質量関数を$${1/13}$$づつ均等に混合したものになります。

※途中の６つの関数を記載省略

■分析・分析結果
分析の仕方や分析数値はテキストの記述が正確だと思いますので、テキストの読み込みをおすすめします。
PyMCの自己流モデルの推論値を利用した分析は「PyMC実装」章をご覧ください。

PyMC実装

Let's enjoy PyMC & Python !

準備

１．インポート

# インポート

# 数値・確率計算
import pandas as pd
import numpy as np

# PyMC
import pymc as pm
import arviz as az

# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'

# ワーニング表示の抑制
import warnings
warnings.simplefilter('ignore')

Sのデータ準備

性行為 S の経験人数に関する回答値を確認します。

１．データの読み込み
テキストの R スクリプト内に記述された回答値データを設定します。

# データの読み込み ※Rスクリプト掲載データを引用（昇順ソートしています）

data = [-9, -9, -9, -8, -8, -8, -8, -8, -8, -8, -7, -7, -7, -7, -6, -6, -6,
        -5, -5, -5, -5, -5, -4, -4, -4, -3, -3, -3, -3, -3, -2, -2, -2, -2, -2,
        -1, -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
        4, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 5, 5, 5,
        6, 6, 6, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
        8, 8, 8, 8, 9, 9, 9, 9, 9, 9, 9, 9, 9, 10, 10, 10, 12, 15, 15]

print('標本サイズ: ', len(data))

【実行結果】
標本サイズは130です。

２．データの外観の確認
まずは要約統計量を確認します。

# 要約統計量の表示
pd.DataFrame({'回答値（マスク済み）': data}).describe().round(2)

【実行結果】
携帯番号下１桁（０から９）でマスクされた回答値です。
最小値の$${-9}$$はマスクの最大値９がマイナスされたものと思われるので、おそらく経験人数は０人でしょう。
最大値の$${+15}$$については、マスクが＋９だと経験人数は６人、マスクが－９だと経験人数は２４人。

箱ひげ図で可視化してみましょう。
テキストの図4.1に相当します。

# 箱ひげ図の描画 ※図4.1に相当

# 描画領域の設定
plt.figure(figsize=(3, 4))
# 箱ひげ図の描画
sns.boxplot(data, color='tab:blue', fill=False)
# スウォームプロットの描画
sns.swarmplot(data, color='blue', alpha=0.3)
plt.ylabel('回答値：経験人数＋マスク');

【実行結果】
スウォームプロットを重ねてデータ点の分布を見通せるようにしました。
－９から９くらいまで満遍なく分布しているような印象です。
回答値＝７が最多回答です。

テキストの表4.1に示された回答値$${y}$$の平均値と標準誤差を計算してみます。

# 平均値と標準誤差の計算 「不偏分散 ÷ 標本サイズ」の平方根 ※表4.1のyに相当
mean1 = np.mean(data)
se1 = np.sqrt(np.var(data, ddof=1) / len(data))
print(f'平均:{mean1:.3f}, 標準誤差:{se1:.3f}')

【実行結果】

Sのポアソン分布モデル構築

モデルの数式表現
目指したいPyMCのモデルの雰囲気を混ぜた「なんちゃって数式」表記です。

$$
\begin{align*}
\lambda &\sim \text{Uniform}\ (\text{lower}=0,\ \text{upper}=100) \\
components &= [ \\
&\text{Poisson.dist}\ (\text{mu}=\lambda)-9, \\
&\text{Poisson.dist}\ (\text{mu}=\lambda)-8, \\
&\quad \vdots \\
&\text{Poisson.dist}\ (\text{mu}=\lambda)-0, \\
&\text{Poisson.dist}\ (\text{mu}=\lambda)+0, \\
&\quad \vdots \\
&\text{Poisson.dist}\ (\text{mu}=\lambda)+8, \\
&\text{Poisson.dist}\ (\text{mu}=\lambda)+9, \\
&] \\
likelihood &\sim \text{Mixture}\ (\text{w}=w,\ \text{comp\_dist}=components)\\
\end{align*}
$$

１．初期値設定
混合分布の混合比率$${w}$$を算出します。
回答値$${y}$$の値に応じてマスクの数$${N}$$を求めます。
マスクの取りうる値が$${-9}$$から$${9}$$まで20個あるので、重み$${w}$$の配列は１行あたり20の要素で構成しています。
コードでは、for文で配列の先頭から$${N}$$個までに$${1/N}$$をセットします。
残りの要素は$${0}$$とし、重み０の混合分布が使われないようにします。

# 初期値設定

# データのインデックス、カテゴリ変数のインデックスの取得
idx_data = list(range(len(data)))

# 混合分布の混合比率w：観測データ（全130行）ごとにマスク個数20個分の均等割合を算出
# wの初期化
w = np.zeros((len(data), 20))
# 観測値yごとに可能なマスク個数を算出
num_masks = [y+10 if y<=-1 else y+11 if y<=9 else 20 for y in data]
# w(numpy配列)の作成
for i, num_mask in enumerate(num_masks):
    # 配列先頭からマスクの個数までの要素を1/Nにする
    w[i, :num_mask] = 1 / num_mask
# wの表示
print('混合比率w: shape=', w.shape)
print(w)

【実行結果】
混合比率$${w}$$の一部を表示しています。

２．モデルの定義
数式表現を実直にモデル記述します。
ポアソン分布＋マスクｍは、カスタム分布クラス pm.CustomDist.dist クラス で実装しています。

# モデルの定義

# logp（対数尤度）関数の定義 ★変更点
def custom_logp1(value, lam, m):
    '''
    value: 観測値 (PyTensor変数が渡される)
    lam: パラメータ
    m: マスク
    '''
    # 「観測値 value から m を引いた値がポアソン分布(λ)に従う」という設定
    return pm.logp(pm.Poisson.dist(mu=lam), value - m)

# coordsの定義
coords = {'data': idx_data}

# モデルの定義
with pm.Model(coords=coords) as model1:

    # Dataの定義：目的変数
    y = pm.Data('y', value=data, dims='data')
    
    # 事前分布の定義
    lam = pm.Uniform('lam', lower=0, upper=100)
    
    # 混合分布の構成要素の定義
    # マスクのfull集合 [-9, -8, …, 0, 0, …, 9]
    masks = np.hstack([np.arange(-9, 1), np.arange(0, 10)])
    # リスト内包表記で20個(マスク個数分)のカスタムポアソン分布を生成 ★変更点
    components = [
        pm.CustomDist.dist(
            lam, m, 
            logp=custom_logp1,            # logp関数を指定
            dist=lambda lam, m, size: (   # 標本サンプリング関数を指定
                pm.Poisson.dist(mu=lam, size=size) + m
            )
        ) for m in masks
    ]
    
    # 尤度：混合分布
    likelihood = pm.Mixture('likelihood', w=w, comp_dists=components,
                            observed=y, dims='data'))

【モデル注釈】

• coordの定義
  座標に名前を付けたり、その座標が取りうる値を設定できます。
  今回は次の１つを設定しました。

  • データの行の座標：名前「data」、値「行インデックス」

• dataの定義
  回答値$${y}$$を設定しました。

• パラメータの事前分布

  • キー項目の内容に関する経験人数$${\lambda}$$の事前分布は区間$${[0,100]}$$の一様分布です。

• 尤度

  • 20のポアソン分布の混合分布です。

  • 混合分布は次のように実装しました。

    • 尤度likelihoodで 混合分布「pm.Mixture」を指定。
      混合分布の構成は「components」を参照。
      各分布の重み（混合比率）はwで指定。

    • componentsのリストでは以下の設定を用いて、リスト内包表記で20個のカスタム分布の確率分布「pm.Custom.dist()」を生成。

      • 対数確率 logp にて、「回答値 $${y}$$ とマスクｍの差」に対するポアソン分布およびパラメータ $${\lambda}$$ を指定。

      • 分布 dist にて Python の無名関数 lambda を利用して、「ポアソン分布＋マスクｍ」による回答値 $${y}$$ の標本サンプリング関数を指定。

logp を用いることで、以前の記事で発生したエラーを回避できました！
そして logp を手に入れてカスタム分布の表現力がパワーアップしました！

なお、パラメータ $${\lambda}$$ の推定値は、「回答値 $${y}$$ とマスクｍの差」、つまり、「Sの経験人数」の平均値です。

３．モデルの外観の確認

# モデルの表示
model1

【実行結果】
２行目の混合分布は右が見切れています。
20個のカスタム分布（ポアソン分布＋マスク）で構成しています。

# モデルの可視化
pm.model_to_graphviz(model1)

【実行結果】
図はとてもシンプルな構造です。

４．事後分布からのサンプリング
MCMCを実行します。
NUTSサンプラーには「nutpie」を利用しました。
nutpie公式サイト▶ https://pymc-devs.github.io/nutpie/
処理時間はおよそ50秒でした。

# 事後分布からのサンプリング ※NUTSサンプラーにnutpieを使用
# テキスト：iter=21000, warmup=1000, chains=5
with model1:
    idata1 = pm.sample(draws=20000, tune=1000, chains=5,
                       nuts_sampler='nutpie', random_seed=1234)

【実行結果】
以前の記事では、このMCMCサンプリング処理でエラー発生しました…
今回は、エラー発生が無く、Divergences の発生も無く、無事に処理できました！

５．サンプリングデータの確認
$${\hat{R}}$$、事後分布の要約統計量、トレースプロットを確認します。
事後分布の収束確認はテキストにならって $${\hat{R}≤1.1}$$ としています。

# r_hat>1.01の確認

# 設定
idata_in = idata1        # idata名
threshold = 1.01         # しきい値

# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())

【実行結果】
ひとまず、しきい値を $${\hat{R}>1.01}$$ で確認しました。
$${\hat{R}>1.01}$$ のパラメータは「0」件です。
全てのパラメータが $${\hat{R} \leq 1.01}$$ であることを確認できました。

ざっくり事後分布サンプルの要約統計量とトレースプロットを確認します。

# 推論データの要約統計情報の表示
var_names = ['lam']
pm.summary(idata1, hdi_prob=0.95, var_names=var_names, round_to=3)

【実行結果】
Sの経験人数の平均パラメータ $${\lambda}$$ の事後平均は $${1.547}$$ 人です。

# トレースプロットの表示
pm.plot_trace(idata1, var_names=var_names, compact=False, figsize=(8, 2))
plt.tight_layout();

【実行結果】
右のグラフはまんべんなく描画されていて、収束している感じがします。
左のグラフの５つの色＝５本のマルコフ連鎖chainは重なっていて、ほぼ同一の推論値をはじき出しています。

平均パラメータ $${\lambda}$$ の事後分布を可視化しましょう。
事後分布サンプリングデータでヒストグラムを描画します。
seaborn の histplot でKDE曲線付きのグラフを作成します。

# パラメータの事後分布の可視化 ※図4.3に相当

# 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 4))

# λのヒストグラムの描画
sns.histplot(idata1.posterior.lam.data.flatten(), bins=30, kde=True,
             stat='density', color='steelblue', edgecolor='white', alpha=0.5,
             ax=ax)
ax.set(title='λの事後分布', xlabel='λ', ylabel='density')
plt.show()

【実行結果】
ほぼ左右対称の分布です！

Sのポアソン分布モデルの分析

Sの経験人数について、テキストにならって分析を進めます。

１．事後分布
平均パラメータ $${\lambda}$$ の事後統計量を表示します。
テキストの表4.1に相当します。

# 推論データの事後統計量の算出 ※表4.1に相当

# 要約統計量計算関数の定義
def calc_stats(x):
    ci = np.quantile(x, q=[0.025, 0.975])
    hdi = az.hdi(x, hdi_prob=0.95)
    return np.mean(x), np.std(x), ci[0], ci[1], hdi[0], hdi[1]

# 要約統計量データフレームの作成
stats_df1 = pd.DataFrame(
    {'λ': calc_stats(idata1.posterior.lam.data.flatten())},
    index = ['EAP', 'post.sd', 'CI2.5%', 'CI97.5%', 'HDI2.5%', 'HDI97.5%']
    ).T

# データフレームの表示
display(stats_df1.round(3))

【実行結果】
テキストとほぼ同じ結果になりました！

【分析】
パラメータ $${\lambda}$$ の平均値（EAP）は $${1.547}$$、事後標準偏差（post.sd）は $${0.297}$$、95%信用区間は $${[1.008,\ 2.171]}$$ です。
事後分布に基づき、経験人数は 95% の確信で $${1.008}$$ 人から $${2.171}$$ 人である、と言えます。

２．事後予測分布
続いて、事後予測分布を見ましょう。
今回は、MCMCサンプルを用いて、numpy の ポアソン乱数による 経験人数 $${x^*}$$ の事後予測サンプルを生成します。
$${x^*}$$ はポアソン分布の不確実性（ばらつき）を含んでいます。

図 4.4 に相当する事後予測プロットを描画します。

# 事後予測分布の可視化 ※図4.4に相当

# 推定された lam のサンプリング結果（全サンプル）を取得
post_lam = idata1.posterior['lam'].values.flatten()

# 推定された lam を使って、マスクなしの純粋なポアソン乱数を生成
rng = np.random.default_rng(seed=1234)  # 乱数生成器の初期化
pure_y_samples = rng.poisson(post_lam)

# 事後予測分布のヒストグラムの描画
plt.hist(pure_y_samples, bins=range(0, max(pure_y_samples) + 3),
         rwidth=0.3, align='left', color='tab:blue')
plt.title('Sの経験人数の事後予測分布')
plt.xlabel('経験人数の事後予測値', fontsize=12)
plt.ylabel('頻度', fontsize=12)
plt.show()

【実行結果】
経験人数は $${0}$$ ～ $${2}$$ 人が多いことが分かります！

表 4.2 の事後予測分布と累積分布を算出しましょう。

# 事後予測分布と累積分布 ※表4.2に相当
unique, counts = np.unique(pure_y_samples, return_counts=True)
ratios = counts / counts.sum()
pd.DataFrame(
    {'x*': ratios, '累積分布': np.cumsum(ratios)}, index=unique
).T.round(2)

【実行結果】

$${x^*}$$ の平均と不偏分散を確認します。

# x*の平均と不偏分散の算出
print(f'x*の平均　　: {np.mean(pure_y_samples):.3f}')
print(f'x*の不偏分散: {np.var(pure_y_samples, ddof=1):.3f}')

【実行結果】
テキストと比べて不偏分散がやや大きい感じがします。
（テキストの不偏分散は $${1.278}$$）

３．PHC曲線：研究仮説が正しい確率
最後に、事後分布のMCMCサンプルを用いてPHC 曲線を描きます。
PHC曲線は、横軸にしきい値 $${c}$$ を置き、$${c}$$ よりも大きいMCMCサンプルの比率をプロットするものです。
平均的な経験人数 $${\lambda}$$ が $${c}$$ 人よりも大きい確率を示します。

PHC曲線を描画する関数を定義します。

# 推論データidataを受け取ってPHC曲線を描画する関数

def phc_plot(data, xrange=10):
    
    ## データ加工
    # x軸目盛値となるcの値を生成
    x = np.linspace(0, xrange, 1001)
    # λ, μ > c の割合を計算
    y = [(data > c).sum() / len(data) for c in x]
    
    ## 描画処理
    plt.plot(x, y)
    plt.xlabel('c')
    plt.ylabel('確率:PHC')

PHC曲線を描画します。
図 4.5 に相当します。

# PHC曲線の描画 ※図4.5に相当
plt.figure(figsize=(6, 4))
phc_plot(idata1.posterior.lam.data.flatten(), xrange=5)
plt.title(f'「平均値 λ > c 」の確率：PHC曲線')
plt.legend(['経験人数のPHC'])
plt.grid(lw=0.3);

【実行結果】
経験人数が１～２人の間で確率が急降下しています。

テキストにならって、$${c=1}$$：１人よりも多い確率と $${c=2}$$：２人よりも多い確率を計算します。

# Sの経験人数の平均が1より大きい（p>1）というPHCの算出
p = idata1.posterior.lam.data.flatten()
c = 1
print(f'p > {c}の確率: {(p > c).sum() / len(p):.1%}')

【実行結果】
平均経験人数が１人よりも多い確率は $${97.7\%}$$ です。

# Sの経験人数の平均が2より大きい（p>2）というPHCの算出
p = idata1.posterior.lam.data.flatten()
c = 2
print(f'p > {c}の確率: {(p > c).sum() / len(p):.1%}')

【実行結果】
平均経験人数が２人よりも多い確率は $${7.1\%}$$ です。

ビジュアルで「経験人数は少なくとも１人はいるだろう」と直感的な理解を促進できます！
（平均経験人数の事後平均 EAP は $${1.547}$$ 人でした。）

S の分析は以上です。

Cのデータ準備

続いて、カンニング C の経験回数のモデリングに進みます。
テキストにならって「ポアソン分布モデル」と「正規分布モデル」に取り組みます。

なお、Cの経験回数のAR法では、マスクを「誕生月」としています。
マスクの要素は$${-12, -11, \cdots, -1, 1, \cdots, 11, 12}$$となり、要素数は24です。

１．データの読み込み
テキストの R スクリプト内に記述された回答値データを設定します。

# データの読み込み ※Rスクリプト掲載データを引用（昇順ソートしています）

data2 = [-11, -8, -7, -7, -7, -7, -7, -7, -7, -6, -5, -5, -5, -5, -5,
         -4, -4, -4, -4, -4, -4, -3, -3, -3, -3, -3, -3, -3,
         -2, -2, -2, -2, -2, -1, -1, -1, -1, -1, -1,
         0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1,
         2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3, 3,
         4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 6, 6,
         7, 7, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 9, 9, 9, 10, 10, 10, 10,
         11, 12, 16, 18, 19, 19, 20, 24]

print('標本サイズ: ', len(data2))

【実行結果】
標本サイズは125です。

２．データの外観の確認
まずは要約統計量を確認します。

# 要約統計量の表示
pd.DataFrame({'回答値（マスク済み）': data2}).describe().round(2)

【実行結果】
誕生月（１から１２）でマスクされた回答値です。
最小値の$${-11}$$はマスクのうち１１または１２がマイナスされたものと思われるので、おそらく経験回数は０回または１回でしょう。
最大値の$${+24}$$については、マスクが＋１２だと経験回数は１２回、マスクが－１２だと経験回数は３６回。

箱ひげ図で可視化してみましょう。
テキストの図4.6に相当します。

# 箱ひげ図の描画 ※図4.6に相当

# 描画領域の設定
plt.figure(figsize=(3, 4))
# 箱ひげ図の描画
sns.boxplot(data2, color='tab:blue', fill=False)
# スウォームプロットの描画
sns.swarmplot(data2, color='blue', alpha=0.3)
plt.ylabel('回答値：経験回数＋マスク');

【実行結果】
スウォームプロットを重ねてデータ点の分布を見通せるようにしました。
－８から１２くらいまでに概ね収まっていて、ひげの外にはみ出た外れ値的な値も見られます。
回答値＝０が最多回答です。
回答値が０になるケースは、マスクが$${-1}$$から$${-12}$$なので、経験回数が$${1}$$から$${12}$$の場合が該当します。

テキストの表4.3に示された回答値$${y}$$の平均値と標準誤差を計算してみます。

# 平均値と標準誤差の計算 「不偏分散 ÷ 標本サイズ」の平方根 ※表4.3のyに相当
mean2 = np.mean(data2)
se2 = np.sqrt(np.var(data2, ddof=1) / len(data2))
print(f'平均:{mean2:.3f}, 標準誤差:{se2:.3f}')

【実行結果】
平均$${2.272}$$は従来法による経験回数の推定値です。

Cのモデル検討

続いて、Sのモデルを転用してC用にアレンジします。

モデルの数式表現
目指したいPyMCのモデルの雰囲気を混ぜた「なんちゃって数式」表記です。
ポアソン分布に代えて正規分布を用いています。

$$
\begin{align*}
\mu &\sim \text{Uniform}\ (\text{lower}=0,\ \text{upper}=100) \\
\sigma &\sim \text{Uniform}\ (\text{lower}=0,\ \text{upper}=100) \\
components &= [ \\
&\text{Normal.dist}\ (\text{mu}=\mu,\ \text{sigma}=\sigma)-12, \\
&\text{Normal.dist}\ (\text{mu}=\mu,\ \text{sigma}=\sigma)-11, \\
&\quad \vdots \\
&\text{Normal.dist}\ (\text{mu}=\mu,\ \text{sigma}=\sigma)-1, \\
&\text{Normal.dist}\ (\text{mu}=\mu,\ \text{sigma}=\sigma)+1, \\
&\quad \vdots \\
&\text{Normal.dist}\ (\text{mu}=\mu,\ \text{sigma}=\sigma)+11, \\
&\text{Normal.dist}\ (\text{mu}=\mu,\ \text{sigma}=\sigma)+12, \\
&] \\
likelihood &\sim \text{Mixture}\ (\text{w}=w,\ \text{comp\_dist}=components)\\
\end{align*}
$$

初期値設定
混合分布の混合比率$${w}$$を算出します。
回答値$${y}$$の値に応じてマスクの数$${N}$$を求めます。
マスクの取りうる値が$${-12}$$から$${12}$$まで24個あるので、重み$${w}$$の配列は１行あたり24の要素で構成しています。
コードでは、for文で配列の先頭から$${N}$$個までに$${1/N}$$をセットします。
残りの要素は$${0}$$とし、重み０の混合分布が使われないようにします。

# 初期値設定

# データのインデックス、カテゴリ変数のインデックスの取得
idx_data = list(range(len(data2)))

# 混合分布の混合比率w：観測データ（全125行）ごとにマスク個数24個分の均等割合を算出
# wの初期化
w = np.zeros((len(data2), 24))
# 観測値yごとに可能なマスク個数を算出
num_masks = [y+13 if y<=-1 else y+12 if y<=12 else 24 for y in data2]
# w(numpy配列)の作成
for i, num_mask in enumerate(num_masks):
    # 配列先頭からマスクの個数までの要素を1/Nにする
    w[i, :num_mask] = 1 / num_mask
# wの表示
print('混合比率w: shape=', w.shape)
print(w)

【実行結果】
混合比率$${w}$$の一部を表示しています。

Cのポアソン分布モデルの構築

最初にポアソン分布モデルのモデリングを行います。
さきほどの S のポアソン分布モデルのコードを流用して作成します。

１．モデルの定義

# モデルの定義 ポアソン分布モデル

# logp（対数尤度）関数の定義
def custom_logp2(value, lam, m):
    '''
    value: 観測値 (PyTensor変数が渡される)
    lam: パラメータ
    m: マスク
    '''
    # 「観測値 value から m を引いた値がポアソン分布(λ)に従う」という設定
    return pm.logp(pm.Poisson.dist(mu=lam), value - m)

# coordsの定義
coords = {'data': idx_data}

# モデルの定義
with pm.Model(coords=coords) as model2:

    # Dataの定義：目的変数
    y = pm.Data('y', value=data2, dims='data')
    
    # 事前分布の定義
    lam = pm.Uniform('lam', lower=0, upper=100)
    
    # 混合分布の構成要素の定義
    # マスクのfull集合 [-12, -1, …, -1, 1, …, 12]
    masks = np.delete(np.arange(-12, 13), 12)
    # リスト内包表記で24個(マスク個数分)のカスタムポアソン分布を生成
    components = [
        pm.CustomDist.dist(
            lam, m, 
            logp=custom_logp2,            # logp関数を指定
            dist=lambda lam, m, size: (  # 標本サンプリング関数を指定
                pm.Poisson.dist(mu=lam, size=size) + m
            )
        ) for m in masks
    ]
    
    # 尤度：混合分布
    likelihood = pm.Mixture('likelihood', w=w, comp_dists=components,
                            observed=y, dims='data')

２．モデルの外観の確認

# モデルの表示
model2

【実行結果】
２行目の混合分布は右が見切れています。
24個のカスタム分布（ポアソン分布＋マスク）で構成しています。

# モデルの可視化
pm.model_to_graphviz(model2)

【実行結果】
先ほどの S のポアソン分布と同じ顔付きです。

３．事後分布からのサンプリング
処理時間はおよそ１分でした。

# 事後分布からのサンプリング ※NUTSサンプラーにnutpieを使用
# テキスト：iter=21000, warmup=1000, chains=5
with model2:
    idata2 = pm.sample(draws=20000, tune=1000, chains=5,
                       nuts_sampler='nutpie', random_seed=1234)

【実行結果】
Divergences は 0 個です。

４．サンプリングデータの確認
$${\hat{R}}$$、事後分布の要約統計量、トレースプロットを確認します。
事後分布の収束確認はテキストにならって$${\hat{R} \leq 1.1}$$とします。

# r_hat>1.01の確認

# 設定
idata_in = idata2        # idata名
threshold = 1.01         # しきい値

# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())

【実行結果】
ひとまず、しきい値を$${\hat{R} > 1.01}$$で確認しました。
$${\hat{R} > 1.01}$$のパラメータは「0」件です。
全てのパラメータが$${\hat{R} \leq 1.01}$$であることを確認できました。

ざっくり事後分布サンプリングデータの要約統計量とトレースプロットを確認します。

# 推論データの要約統計情報の表示
var_names = ['lam']
pm.summary(idata2, hdi_prob=0.95, var_names=var_names, round_to=3)

【実行結果】
Cの経験回数の平均パラメータ $${\lambda}$$ の事後平均（EAP）は5.008です。

# トレースプロットの表示
pm.plot_trace(idata2, var_names=var_names, compact=False, figsize=(8, 2))
plt.tight_layout();

【実行結果】
右のグラフはまんべんなく描画されていて、収束している感じがします。
左のグラフの５つの色＝５本のマルコフ連鎖chainは重なっていて、ほぼ同一の推論値をはじき出しています。

Cのポアソン分布モデルの分析

１．事後分布
最初にテキストの表 4.3 に相当する事後分布の要約統計量を眺めます。

# 推論データの事後統計量の算出 ※表4.3に相当

# 要約統計量データフレームの作成
stats_df2 = pd.DataFrame(
    {'λ': calc_stats(idata2.posterior.lam.data.flatten())},
    index = ['EAP', 'post.sd', 'CI2.5%', 'CI97.5%', 'HDI2.5%', 'HDI97.5%']
    ).T

# データフレームの表示
display(stats_df2.round(3))

【実行結果】

【分析】
パラメータ $${\lambda}$$ の平均値（EAP）は $${5.008}$$、事後標準偏差（post.sd）は $${0.537}$$、95%信用区間は $${[3.972,\ 6.079]}$$ です。
事後分布に基づき、経験回数は 95% の確信で $${3.972}$$ 回から $${6.079}$$ 回である、と言えます。

事後分布を可視化しましょう。

# パラメータの事後分布の可視化

# 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 4))

# λのヒストグラムの描画
sns.histplot(idata2.posterior.lam.data.flatten(), bins=30, kde=True,
             stat='density', color='steelblue', edgecolor='white', alpha=0.5,
             ax=ax)
ax.set(title='λの事後分布', xlabel='λ', ylabel='density')
plt.show()

【実行結果】
美しいベル型です！

２．事後予測分布
MCMCサンプルを用いて、経験回数 $${x^*}$$ の事後予測分布のサンプル生成と描画をしましょう。

# 事後予測分布の可視化

# 推定された lam のサンプリング結果（全サンプル）を取得
post_lam2 = idata2.posterior['lam'].values.flatten()

# 推定された lam を使って、マスクなしの純粋なポアソン乱数を生成
rng = np.random.default_rng(seed=1234)  # 乱数生成器の初期化
pure_y_samples2 = rng.poisson(post_lam2)

# 事後予測分布のヒストグラムの描画
plt.hist(pure_y_samples2, bins=range(0, max(pure_y_samples2) + 2),
         rwidth=0.3, align='left', color='tab:blue')
plt.title('Cの経験人数の事後予測分布')
plt.xlabel('経験人数の事後予測値', fontsize=12)
plt.ylabel('頻度', fontsize=12)
plt.show()

【実行結果】
４～５回にピークがあります。

$${x^*}$$ の事後予測の平均と不偏分散を確認します。

# 事後予測分布と累積分布
unique, counts = np.unique(pure_y_samples2, return_counts=True)
ratios = counts / counts.sum()
pd.DataFrame(
    {'x*': ratios, '累積分布': np.cumsum(ratios)}, index=unique
).T.round(2)

【実行結果】

３．PHC曲線
テキストには無いですが、ポアソン分布モデルのPHC曲線を見ておきます。

# PHC曲線の描画
plt.figure(figsize=(6, 4))
phc_plot(idata2.posterior.lam.data.flatten(), xrange=10)
plt.title(f'「平均値 λ > c 」の確率：PHC曲線')
plt.legend(['経験回数のPHC'])
plt.grid(lw=0.3);

【実行結果】
経験回数が４～５回の間で確率が急降下しています。

$${c=4}$$：４回よりも多い確率と $${c=6}$$：６回よりも多い確率を計算します。

# Cの経験回数の平均が4より大きい（p>4）というPHCの算出
p = idata2.posterior.lam.data.flatten()
c = 4
print(f'p > {c}の確率: {(p > c).sum() / len(p):.1%}')

【実行結果】
平均経験回数が４回よりも多い確率は $${97.2\%}$$ です。

# Cの経験回数の平均が6より大きい（p>6）というPHCの算出
p = idata2.posterior.lam.data.flatten()
c = 6
print(f'p > {c}の確率: {(p > c).sum() / len(p):.1%}')

【実行結果】
平均経験回数が６人よりも多い確率は $${3.4\%}$$ です。

ビジュアルで経験回数が少なくとも４回はありそうだ、と推測できますね！
（平均経験回数の事後平均 EAP は $${5.008}$$ 回でした。）

Cの正規分布モデルの構築

次に正規分布モデルのモデリングを行います。
ポアソン分布モデルのコードのパターンを流用して、正規分布用に書き換えます。

１．モデルの定義

# モデルの定義 正規分布モデル

# logp（対数尤度）関数の定義
def custom_logp3(value, mu, sigma, m):
    '''
    value: 観測値 (PyTensor変数が渡される)
    mu: 平均
    sigma: 標準偏差
    m: マスク
    '''
    # 「観測値 value から m を引いた値が正規分布(μ, σ)に従う」という設定
    return pm.logp(pm.Normal.dist(mu=mu, sigma=sigma), value - m)

# coordsの定義
coords = {'data': idx_data}

# モデルの定義
with pm.Model(coords=coords) as model3:
    
    # Dataの定義
    y = pm.Data('y', value=data2, dims='data')
  
    # 事前分布の定義
    mu = pm.Uniform('mu', lower=0, upper=100)
    sigma = pm.Uniform('sigma', lower=0, upper=100)
    
    # 混合分布の構成要素の定義
    # マスクのfull集合 [-12, -1, …, -1, 1, …, 12]
    masks = np.delete(np.arange(-12, 13), 12)
    # リスト内包表記で24個(マスク個数分)のカスタム正規分布を生成
    components = [
        pm.CustomDist.dist(
            mu, sigma, m,
            logp=custom_logp3,                   # logp関数を指定
            dist=lambda mu, sigma, m, size: (    # 標本サンプリング関数を指定
                pm.Normal.dist(mu=mu, sigma=sigma, size=size) + m
            )
        ) for m in masks
    ]
    
    # 尤度：混合分布
    likelihood = pm.Mixture('likelihood', w=w, comp_dists=components,
                             observed=y, dims='data')

【モデル注釈】

• coordの定義
  座標に名前を付けたり、その座標が取りうる値を設定できます。
  今回は次の１つを設定しました。

  • データの行の座標：名前「data」、値「行インデックス」

• dataの定義
  回答値{y}$$を設定しました。

• パラメータの事前分布

  • キー項目の内容に関する経験回数の平均$${\mu}$$と標準偏差$${\sigma}$$の事前分布は、区間$${[0,100]}$$の一様分布です。

• 尤度

  • 24の正規分布の混合分布です。

  • 混合分布は次のように実装しました。

    • 尤度likelihoodで 混合分布「pm.Mixture」を指定。
      混合分布の構成は「components」を参照。
      各分布の重み（混合比率）はwで指定。

    • componentsのリストでは以下の設定を用いて、リスト内包表記で24個のカスタム分布の確率分布「pm.Custom.dist()」を生成。

      • 対数確率 logp にて、「回答値 $${y}$$ とマスクｍの差」に対する正規分布およびパラメータ $${\mu,\ \sigma}$$ を指定。

      • 分布 dist にて Python の無名関数 lambda を利用して、「正規分布＋マスクｍ」による回答値 $${y}$$ の標本サンプリング関数を指定。

２．モデルの外観の確認

# モデルの表示
model3

【実行結果】
２行目の混合分布は右が見切れています。
24個のカスタム分布（正規分布＋マスク）で構成しています。

# モデルの可視化
pm.model_to_graphviz(model3)

【実行結果】
隠れミッキーがいます！？

３．事後分布からのサンプリング
処理時間はおよそ１分でした。

# 事後分布からのサンプリング ※NUTSサンプラーにnutpieを使用
# テキスト：iter=21000, warmup=1000, chains=5
with model3:
    idata3 = pm.sample(draws=20000, tune=1000, chains=5,
                       nuts_sampler='nutpie', random_seed=1234)

【実行結果】
Divergences は 0 個です。

４．サンプリングデータの確認
$${\hat{R}}$$、事後分布の要約統計量、トレースプロットを確認します。
事後分布の収束確認はテキストにならって$${\hat{R} \leq 1.1}$$とします。

# r_hat>1.01の確認
# 設定
idata_in = idata3        # idata名
threshold = 1.01         # しきい値

# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())

【実行結果】
ひとまず、しきい値を$${\hat{R} > 1.01}$$で確認しました。
$${\hat{R} > 1.01}$$のパラメータは「0」件です。
全てのパラメータが$${\hat{R} \leq 1.01}$$であることを確認できました。

ざっくり事後分布サンプリングデータの要約統計量とトレースプロットを確認します。

# 推論データの要約統計情報の表示
var_names = ['mu', 'sigma']
pm.summary(idata3, hdi_prob=0.95, var_names=var_names, round_to=3)

【実行結果】
Cの経験回数の平均$${\mu}$$の平均値（テキストのEAP）は5.340です。

# トレースプロットの表示
pm.plot_trace(idata3, var_names=var_names, compact=False, figsize=(8, 4))
plt.tight_layout();

【実行結果】
右のグラフはまんべんなく描画されていて、収束している感じがします。
左のグラフの５つの色＝５本のマルコフ連鎖chainは重なっていて、ほぼ同一の推論値をはじき出しています。

Cの正規分布モデルの分析

Cの経験回数について、正規分布モデルの事後分布に基づいて分析を進めます。

１．事後分布
平均$${\mu}$$と標準偏差$${\sigma}$$の事後統計量を表示します。
テキストの表4.4に相当します。

# 推論データの事後統計量の算出 ※表4.4に相当

# 要約統計量データフレームの作成
stats_df2 = pd.DataFrame(
    {'μ': calc_stats(idata3.posterior.mu.data.flatten()),
     'σ': calc_stats(idata3.posterior.sigma.data.flatten())},
    index = ['EAP', 'post.sd', 'CI2.5%', 'CI97.5%', 'HDI2.5%', 'HDI97.5%']
    ).T

# データフレームの表示
display(stats_df2.round(3))

【実行結果】
テキストとほぼ同じ結果になりました！

【分析】
平均パラメータ$${\mu}$$の平均値（EAP）は$${5.340}$$、事後標準偏差（post.sd）は$${0.516}$$です。

点推定・区間推定を比較します。

$$
\begin{array}{c|ccc}
モデル & 平均値 & 標準偏差・誤差 & 95\% \text{CI}下側 & 95\% \text{CI}上側 \\
\hline
正規分布 & 5.340 & 0.516 & 4.353 & 6.376 \\
ポアソン分布 & 5.008 & 0.537 & 3.972 & 6.079 \\
従来法 & 2.272 & 0.551 & - & - \\
\end{array}
$$

• 平均値に関して

  • 正規分布モデルとポアソン分布モデルは同程度です。

  • 一方で、２つのモデルと従来法の推定値はかけ離れています。

  • 従来法の平均値$${2.272}$$は、正規分布モデル・ポアソン分布モデルの95%CIの区間に含まれていません。

• 標準偏差・誤差・95%信用区間に関して

  • 正規分布モデル・ポアソン分布モデルの標準偏差は、従来法の標準誤差$${0.551}$$よりも小さくなっています。

  • 正規分布モデルの標準偏差、95%信用区間の幅は、ポアソン分布モデルよりも小さくて狭いです。

正規分布モデルの事後分布を可視化します。

# パラメータの事後分布の可視化

# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4), sharey=True)
# μのヒストグラムの描画
sns.histplot(idata3.posterior.mu.data.flatten(), bins=30, kde=True,
             stat='density', color='steelblue', edgecolor='white', alpha=0.5,
             ax=ax1)
ax1.set(title='μの事後分布', xlabel='μ', ylabel='density')
# σのヒストグラムの描画
sns.histplot(idata3.posterior.sigma.data.flatten(), bins=30, kde=True,
             stat='density', color='tomato', edgecolor='white', alpha=0.5,
             ax=ax2)
ax2.set(title='σの事後分布', xlabel='σ')
# 最終化
plt.tight_layout()
plt.show()

【実行結果】
美しいベル型です！

２．事後予測分布
MCMCサンプルを用いて、経験回数 $${x^*}$$ の事後予測分布のサンプル生成と描画をしましょう。

# 事後予測分布の可視化

# 推定された mu, sigma のサンプリング結果（全サンプル）を取得
post_mu = idata3.posterior['mu'].values.flatten()
post_sigma = idata3.posterior['sigma'].values.flatten()

# 推定された mu, sigma を使って、マスクなしの純粋な正規乱数を生成
rng = np.random.default_rng(seed=1234)  # 乱数生成器の初期化
pure_y_samples3 = rng.normal(loc=post_mu, scale=post_sigma)

# 事後予測分布のヒストグラムの描画
plt.hist(pure_y_samples3, bins=30, align='left', color='tab:blue',
         edgecolor='white')
plt.title('Cの経験人数の事後予測分布')
plt.xlabel('経験人数の事後予測値', fontsize=12)
plt.ylabel('頻度', fontsize=12)
plt.show()

【実行結果】
４～６回にピークがあります。
負の値が発生していて、ちょっと残念な気持ちです…

$${x^*}$$ の事後予測の平均と不偏分散を確認します。

# x*の平均と不偏分散の算出
print(f'x*の平均　　: {np.mean(pure_y_samples3):.3f}')
print(f'x*の不偏分散: {np.var(pure_y_samples3, ddof=1):.3f}')

【実行結果】

３．PHC曲線
続いてPHC曲線を描画します。
テキストの図4.7に相当し、正規分布モデルを仮定した上での「Cの経験回数の平均が$${c}$$よりも大きい」という研究仮説が正しい確率の曲線です。

# PHC曲線の描画 ※図4.7に相当
plt.figure(figsize=(6, 4))
phc_plot(idata3.posterior.mu.data.flatten())
plt.title(f'「平均値 μ > c 」の確率：PHC曲線')
plt.legend(['経験回数のPHC'])
plt.grid(lw=0.3);

【実行結果】
テキストとほぼ同じ結果を得られました。
４回以上６回以下の間で確率が大きく変動（急降下）していることが分かります。

経験回数の平均が 4 より大きい（$${p>4}$$）というPHC、6 より大きい（$${p>6}$$）というPHCを計算してみます。

# Cの経験回数の平均が4より大きい（p>4）というPHCの算出
p = idata3.posterior.mu.data.flatten()
c = 4
print(f'p > {c}の確率: {(p > c).sum() / len(p):.1%}')

【実行結果】
４回を上回る確率はほぼ100%です。
Cの経験の集団平均は、少なくとも４回はありそうです。

# Cの経験回数の平均が6より大きい（p>6）というPHCの算出
p = idata3.posterior.mu.data.flatten()
c = 6
print(f'p > {c}の確率: {(p > c).sum() / len(p):.1%}')

【実行結果】
６回を上回る確率は小さくなりましたが、$${10\%}$$ あるようです。

idata のファイル保存は省略します。

以上で第４章は終了です。

おわりに

３つの間接質問法記事の起源

第２～４章で３つの間接質問法の手法に取り組みました。
個人的には第３章のランダム回答法がシンプルで分かりやすかったです。
ただランダム回答法が適用できるのは回答値が２値の場合のようです。
カウントデータや連続値なども扱えるシンプルな方法（アイテムカウント法・AR法以外）は無いものかと、Web検索していたら・・・
次のサイトを発見しました！

https://www.waseda.jp/sem-toyoda-lab/backbook/gform.html

第２～４章の執筆者の顔ぶれも。
当時はベイズを利用していなかったようです。
ご参考まで。

🍀🍀🍀

ライブラリのバージョンアップ

この改造記事は 2026年3月9日に書いています。
実は、3月初めに ArviZ が Ver. 1.0.0 になるという、メジャーアップデート（というか正式バージョンの公開）がありました。
気になるのは、大きな変化があったことです。

• ３つのライブラリに分割

  • arviz-base

  • arviz-stats

  • arviz-plots

• InferenceData 形式が発展的解消され、xarray.DataTree 形式に変更

この変更に伴って、PyMC も変更対応が実施される感じです…

いままで書いた（ブログにアップした）PyMC ＆ ArviZ コード、どうしようかな…
悩み多き年度末です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。