見出し画像

StanとRでベイズ統計モデリングをPyMC Ver.5で写経～第７章「7.9 外れ値」

ネイピア DS

第７章「回帰分析の悩みどころ」

書籍の著者　松浦健太郎 先生

この記事は、テキスト第７章「回帰分析の悩みどころ」の7.9節「外れ値」の PyMC5写経 を取り扱います。
コーシー分布と Student の t 分布の２つのモデルに取り組みます。
外れ値を含めて解析を進めます。

はじめに

StanとRでベイズ統計モデリングの紹介

この記事は書籍「StanとRでベイズ統計モデリング」（共立出版、「テキスト」と呼びます）のベイズモデルを用いて、PyMC Ver.5で「実験的」に写経する翻訳的ドキュメンタリーです。

テキストは、2016年10月に発売され、ベイズモデリングのモデル式とプログラミングに関する丁寧な解説とモデリングの改善ポイントを網羅するチュートリアル「実践解説書」です。もちろん素晴らしいです！
「アヒル本」の愛称で多くのベイジアンに愛されてきた書籍です！

テキストに従ってStanとRで実践する予定でしたが、RのStan環境を整えることができませんでした（泣）
そこでこのシリーズは、テキストのベイズモデルをPyMC Ver.5に書き換えて実践します。

引用表記

この記事は、出典に記載の書籍に掲載された文章及びコードを引用し、適宜、掲載文章とコードを改変して書いています。
【出典】
「StanとRでベイズ統計モデリング」初版第13刷、著者 松浦健太郎、共立出版

記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！

PyMC環境の準備

Anacondaを用いる環境構築とGoogle ColaboratoryでPyMCを動かす方法について、次の記事にまとめています。
「PyMCを動かすまでの準備」章をご覧ください。


7.9 外れ値

インポート

### インポート

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

# 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')

データの読み込み

サンプルコードのデータを読み込みます。

### データの読み込み ◆data-outlier.txt
# Y:目的変数, X:説明変数

data = pd.read_csv('./data/data-outlier.txt')
print('data.shape: ', data.shape)
display(data.head())

【実行結果】

要約統計量と相関係数を算出します。

### 要約統計量の表示
data.describe().round(2)

【実行結果】

### 相関係数の表示
data.corr().round(3)

【実行結果】

散布図で可視化します。
テキスト図7.11に相当します。

### 散布図の描画 ◆図7.11

# 外れ値の色を変えるためのデータ作成
hue = data['Y'].index !=2
# 描画領域の設定
plt.figure(figsize=(5, 4))
ax = plt.subplot()
# 散布図の描画
sns.scatterplot(ax=ax, data=data, x='X', y='Y', hue=hue, palette='Pastel1',
                s=100, legend=None)
# 修飾
ax.set(title='散布図', ylim=(-29, 79), yticks=range(-25, 76, 25),
       xticks=range(0, 11))
ax.grid(lw=0.5)
plt.show()

【実行結果】
$${X=3}$$の点が飛び出しています。

Y の予測に用いる X の値を設定します。

### 予測に用いるXの作成
X_new = np.linspace(0, 11, 101)

モデル式7-8 正規分布

モデルの定義です。

### モデルの定義 ◆モデル式7-8

with pm.Model() as model1:
    
    ### データ関連定義
    ## coordの定義
    model1.add_coord('data', values=data.index, mutable=True)
    model1.add_coord('dataNew', values=range(len(X_new)), mutable=True)
    ## dataの定義
    # 目的変数 Y
    Y = pm.ConstantData('Y', value=data['Y'].values, dims='data')
    # 説明変数 X
    X = pm.ConstantData('X', value=data['X'].values, dims='data')
    # 説明変数 XNew ：予測用
    XNew = pm.ConstantData('XNew', value=X_new, dims='dataNew')

    ### 事前分布
    a = pm.Uniform('a', lower=-100, upper=100)
    b = pm.Uniform('b', lower=-100, upper=100)
    sigma = pm.Uniform('sigma', lower=0, upper=100)

    ### 線形予測子
    mu = pm.Deterministic('mu', a + b * X, dims='data')

    ### 尤度関数
    obs = pm.Normal('obs', mu=mu, sigma=sigma, observed=Y, dims='data')

    ### 計算値
    muNew = a + b * XNew
    yNew = pm.Normal('yNew', mu=muNew, sigma=sigma, dims='dataNew')

モデルの定義内容を見ます。

### モデルの表示
model1

【実行結果】

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

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

MCMCを実行します。

### 事後分布からのサンプリング 15秒
with model1:
    idata1  = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.998,
                        nuts_sampler='numpyro', random_seed=1234)

【実行結果】省略

Pythonで事後分布からのサンプリングデータの確認を行います。
Rhatの確認から。
テキストの収束条件は「chainを3以上にして$${\hat{R}<1.1}$$のとき」です。

### r_hat>1.1の確認
# 設定
idata_in = idata1        # idata名
threshold = 1.02         # しきい値

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

【実行結果】
収束条件を満たしています。

事後統計量を表示します。

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

【実行結果】

トレースプロットを描画します。

### トレースプロットの表示
var_names = ['a', 'b', 'sigma', 'mu', 'yNew']
pm.plot_trace(idata1, compact=True, var_names=var_names)
plt.tight_layout();

【実行結果】

パラメータの事後統計量の要約を算出します。

### パラメータの要約を確認

# mean,sd,2.5%,25%,50%,75%,97.5%パーセンタイル点をデータフレーム化する関数の定義
def make_stats_df(y):
    probs = [2.5, 25, 50, 75, 97.5]
    columns = ['mean', 'sd'] + [str(s) + '%' for s in probs]
    quantiles = pd.DataFrame(np.percentile(y, probs, axis=0).T, index=y.columns)
    tmp_df = pd.concat([y.mean(axis=0), y.std(axis=0), quantiles], axis=1)
    tmp_df.columns=columns
    return tmp_df

# 要約統計量の算出・表示
vars = ['a', 'b', 'sigma']
param_samples = idata1.posterior[vars].to_dataframe().reset_index(drop=True)
display(make_stats_df(param_samples).round(2))

【実行結果】

$${Y}$$の観測値と予測値のプロットを描画します。
テキスト図7.12左に相当します。

### 散布図の描画 ◆図7.12左

## 描画用データの作成
# MCMCサンプリングデータからyNewを取り出し
y_pred_samples1 = idata1.posterior.yNew.stack(sample=('chain','draw')).data
# y_predの中央値、95%CI、50%CIを算出
y_pred_median1 = np.median(y_pred_samples1, axis=1)
y_pred_95ci1 = np.quantile(y_pred_samples1, q=[0.025, 0.975], axis=1)
y_pred_50ci1 = np.quantile(y_pred_samples1, q=[0.250, 0.750], axis=1)

## 描画処理
# 描画領域の設定
plt.figure(figsize=(5, 4))
ax = plt.subplot()
# 観測値の散布図の描画
sns.scatterplot(ax=ax, data=data, x='X', y='Y', s=80, alpha=0.7, label='観測値')
# 予測値の中央値の折れ線グラフの描画
ax.plot(X_new, y_pred_median1, color='tab:red', label='予測値(中央値)')
# 予測値の95%CIの描画
ax.fill_between(X_new, y_pred_95ci1[0], y_pred_95ci1[1], color='tomato',
                 alpha=0.2, label='95%CI')
# 予測値の50%CIの描画
ax.fill_between(X_new, y_pred_50ci1[0], y_pred_50ci1[1], color='tomato',
                 alpha=0.5, label='50%CI')
# 修飾
ax.set(title='$y$ の観測値と予測値のプロット\n正規分布',
       ylim=(-29, 79), yticks=range(-25, 76, 25), xticks=range(0, 12))
ax.legend(bbox_to_anchor=(1.41, 1))
ax.grid(lw=0.5);

【実行結果】

モデル式7-9 コーシー分布

モデルの定義です。
$${yNew}$$のMCMCサンプルは生成させず、散布図描画の際に予測分布を算出します。
MCMCサンプルデータには大きなノイズが入ってしまうからです（原因不明）。

### モデルの定義 ◆モデル式7-9, model7-9.stan ※yNewはMCMC外で算出

with pm.Model() as model2:
    
    ### データ関連定義
    ## coordの定義
    model2.add_coord('data', values=data.index, mutable=True)
    model2.add_coord('dataNew', values=range(len(X_new)), mutable=True)
    ## dataの定義
    # 目的変数 Y
    Y = pm.ConstantData('Y', value=data['Y'].values, dims='data')
    # 説明変数 X
    X = pm.ConstantData('X', value=data['X'].values, dims='data')

    ### 事前分布
    a = pm.Uniform('a', lower=-100, upper=100)
    b = pm.Uniform('b', lower=-100, upper=100)
    sigma = pm.Uniform('sigma', lower=0, upper=45)

    ### 線形予測子
    mu = pm.Deterministic('mu', a + b * X, dims='data')

    ### 尤度関数
    obs = pm.Cauchy('obs', alpha=mu, beta=sigma, observed=Y, dims='data')

モデルの定義内容を見ます。

### モデルの表示
model2

【実行結果】

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

【実行結果】

MCMCを実行します。

### 事後分布からのサンプリング 10秒
with model2:
    idata2  = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.9,
                        nuts_sampler='numpyro', random_seed=1234)

【実行結果】省略

Pythonで事後分布からのサンプリングデータの確認を行います。
Rhatの確認から。
テキストの収束条件は「chainを3以上にして$${\hat{R}<1.1}$$のとき」です。

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

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

【実行結果】
収束条件を満たしています。

事後統計量を表示します。

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

【実行結果】

トレースプロットを描画します。

### トレースプロットの表示
var_names = ['a', 'b', 'sigma', 'mu']
pm.plot_trace(idata2, compact=True, var_names=var_names)
plt.tight_layout();

【実行結果】

パラメータの事後統計量の要約を算出します。

### パラメータの要約を確認

# mean,sd,2.5%,25%,50%,75%,97.5%パーセンタイル点をデータフレーム化する関数の定義
def make_stats_df(y):
    probs = [2.5, 25, 50, 75, 97.5]
    columns = ['mean', 'sd'] + [str(s) + '%' for s in probs]
    quantiles = pd.DataFrame(np.percentile(y, probs, axis=0).T, index=y.columns)
    tmp_df = pd.concat([y.mean(axis=0), y.std(axis=0), quantiles], axis=1)
    tmp_df.columns=columns
    return tmp_df

# 要約統計量の算出・表示
vars = ['a', 'b', 'sigma']
param_samples = idata2.posterior[vars].to_dataframe().reset_index(drop=True)
display(make_stats_df(param_samples).round(2))

【実行結果】

$${Y}$$の観測値と予測値のプロットを描画します。
テキスト図7.12左に相当します。
冒頭で$${Y}$$の予測値を算出しています。

### 散布図の描画 ◆図7.12右

## 描画用データの作成１：yの予測値の算出
# MCMCサンプリングデータからa,b,sigmaを取り出し
a_samples2 = idata2.posterior.a.stack(sample=('chain','draw')).data
b_samples2 = idata2.posterior.b.stack(sample=('chain','draw')).data
sigma_samples2 = idata2.posterior.sigma.stack(sample=('chain','draw')).data
# scipyのcahchyで乱数生成してyの予測値を算出 #乱数シードを固定すると値が超乱れる
y_pred_samples2 = np.array([
    stats.cauchy.rvs(loc=a_sample + b_sample * X_new, scale=sigma_sample)
    for (a_sample, b_sample, sigma_sample)
    in zip(a_samples2, b_samples2, sigma_samples2)]).T

## 描画用データの作成２：y_predの中央値、95%CI、50%CIを算出
y_pred_median2 = np.median(y_pred_samples2, axis=1)
y_pred_95ci2 = np.quantile(y_pred_samples2, q=[0.025, 0.975], axis=1)
y_pred_50ci2 = np.quantile(y_pred_samples2, q=[0.250, 0.750], axis=1)

## 描画処理
# 描画領域の設定
plt.figure(figsize=(5, 4))
ax = plt.subplot()
# 観測値の散布図の描画
sns.scatterplot(ax=ax, data=data, x='X', y='Y', s=80, alpha=0.7, label='観測値')
# 予測値の中央値の折れ線グラフの描画
ax.plot(X_new, y_pred_median2, color='tab:red', label='予測値(中央値)')
# 予測値の95%CIの描画
ax.fill_between(X_new, y_pred_95ci2[0], y_pred_95ci2[1], color='tomato',
                 alpha=0.2, label='95%CI')
# 予測値の50%CIの描画
ax.fill_between(X_new, y_pred_50ci2[0], y_pred_50ci2[1], color='tomato',
                 alpha=0.5, label='50%CI')
# 修飾
ax.set(title='$y$ の観測値と予測値のプロット\nコーシー分布',
       ylim=(-29, 79), yticks=range(-25, 76, 25), xticks=range(0, 12))
ax.legend(bbox_to_anchor=(1.41, 1))
ax.grid(lw=0.5);

【実行結果】

正規分布モデルとコーシー分布モデルの比較

$${Y}$$の観測値と予測値のプロットを２つ横並びにして描画します。
テキスト図7.12に相当します。

## 描画横並び ◆図7.12

# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 4), sharex=True, sharey=True)

# 観測値の散布図の描画
sns.scatterplot(ax=ax1, data=data, x='X', y='Y', s=80, alpha=0.7)
# 予測値の中央値の折れ線グラフの描画
ax1.plot(X_new, y_pred_median1, color='tab:red')
# 予測値の95%CIの描画
ax1.fill_between(X_new, y_pred_95ci1[0], y_pred_95ci1[1], color='tomato',
                 alpha=0.2)
# 予測値の50%CIの描画
ax1.fill_between(X_new, y_pred_50ci1[0], y_pred_50ci1[1], color='tomato',
                 alpha=0.5)
# 修飾
ax1.set(title='$y$ の観測値と予測値のプロット\n正規分布',
       ylim=(-29, 79), yticks=range(-25, 76, 25), xticks=range(0, 12))
ax1.grid(lw=0.5);


# 観測値の散布図の描画
sns.scatterplot(ax=ax2, data=data, x='X', y='Y', s=80, alpha=0.7, label='観測値')
# 予測値の中央値の折れ線グラフの描画
ax2.plot(X_new, y_pred_median2, color='tab:red', label='予測値(中央値)')
# 予測値の95%CIの描画
ax2.fill_between(X_new, y_pred_95ci2[0], y_pred_95ci2[1], color='tomato',
                 alpha=0.2, label='95%CI')
# 予測値の50%CIの描画
ax2.fill_between(X_new, y_pred_50ci2[0], y_pred_50ci2[1], color='tomato',
                 alpha=0.5, label='50%CI')
# 修飾
ax2.set(title='$y$ の観測値と予測値のプロット\nコーシー分布',
       ylim=(-29, 79), yticks=range(-25, 76, 25), xticks=range(0, 12))
ax2.grid(lw=0.5)
ax2.legend(bbox_to_anchor=(1.5, 1))
plt.tight_layout();

【実行結果】
コーシー分布の予測値（中央値、ベイズ信用区間）が$${X=3}$$の外れ値に引きずられていない様子が分かります。

モデルExtra Student の t 分布

モデルの定義です。
$${yNew}$$のMCMCサンプルは生成させず、散布図描画の際に予測分布を算出します。
MCMCサンプルデータには大きなノイズが入ってしまうからです（原因不明）。

### モデルの定義 ※yNewはMCMC外で算出

with pm.Model() as model3:
    
    ### データ関連定義
    ## coordの定義
    model3.add_coord('data', values=data.index, mutable=True)
    model3.add_coord('dataNew', values=range(len(X_new)), mutable=True)
    ## dataの定義
    # 目的変数 Y
    Y = pm.ConstantData('Y', value=data['Y'].values, dims='data')
    # 説明変数 X
    X = pm.ConstantData('X', value=data['X'].values, dims='data')

    ### 事前分布
    a = pm.Uniform('a', lower=-100, upper=100)
    b = pm.Uniform('b', lower=-100, upper=100)
    sigma = pm.Uniform('sigma', lower=0, upper=100)
    nu = pm.Uniform('nu', lower=0, upper=10)

    ### 線形予測子
    mu = pm.Deterministic('mu', a + b * X, dims='data')

    ### 尤度関数
    obs = pm.StudentT('obs', nu=nu, mu=mu, sigma=sigma, observed=Y,
                      dims='data')

モデルの定義内容を見ます。

### モデルの表示
model3

【実行結果】

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

【実行結果】

MCMCを実行します。

### 事後分布からのサンプリング 10秒
with model3:
    idata3  = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.9,
                        nuts_sampler='numpyro', random_seed=1234)

【実行結果】省略

Pythonで事後分布からのサンプリングデータの確認を行います。
Rhatの確認から。
テキストの収束条件は「chainを3以上にして$${\hat{R}<1.1}$$のとき」です。

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

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

【実行結果】
収束条件を満たしています。

事後統計量を表示します。

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

【実行結果】

トレースプロットを描画します。

### トレースプロットの表示
var_names = ['a', 'b', 'sigma', 'nu', 'mu']
pm.plot_trace(idata3, compact=True, var_names=var_names)
plt.tight_layout();

【実行結果】

パラメータの事後統計量の要約を算出します。

### パラメータの要約を確認

# mean,sd,2.5%,25%,50%,75%,97.5%パーセンタイル点をデータフレーム化する関数の定義
def make_stats_df(y):
    probs = [2.5, 25, 50, 75, 97.5]
    columns = ['mean', 'sd'] + [str(s) + '%' for s in probs]
    quantiles = pd.DataFrame(np.percentile(y, probs, axis=0).T, index=y.columns)
    tmp_df = pd.concat([y.mean(axis=0), y.std(axis=0), quantiles], axis=1)
    tmp_df.columns=columns
    return tmp_df

# 要約統計量の算出・表示
vars = ['a', 'b', 'sigma', 'nu']
param_samples = idata3.posterior[vars].to_dataframe().reset_index(drop=True)
display(make_stats_df(param_samples).round(2))

【実行結果】

$${Y}$$の観測値と予測値のプロットを描画します。
テキスト図7.12左に相当します。
冒頭で$${Y}$$の予測値を算出しています。

### 散布図の描画

## 描画用データの作成１：yの予測値の算出
# MCMCサンプリングデータからa,b,sigmaを取り出し
nu_samples3 = idata3.posterior.nu.stack(sample=('chain','draw')).data
a_samples3 = idata3.posterior.a.stack(sample=('chain','draw')).data
b_samples3 = idata3.posterior.b.stack(sample=('chain','draw')).data
sigma_samples3 = idata3.posterior.sigma.stack(sample=('chain','draw')).data
# scipyのcahchyで乱数生成してyの予測値を算出 #乱数シードを固定すると値が超乱れる
y_pred_samples3 = np.array([
    stats.t.rvs(df=nu_sample, loc=a_sample + b_sample * X_new, scale=sigma_sample)
    for (nu_sample, a_sample, b_sample, sigma_sample)
    in zip(nu_samples3, a_samples3, b_samples3, sigma_samples3)]).T

## 描画用データの作成２：y_predの中央値、95%CI、50%CIを算出
y_pred_median3 = np.median(y_pred_samples3, axis=1)
y_pred_95ci3 = np.quantile(y_pred_samples3, q=[0.025, 0.975], axis=1)
y_pred_50ci3 = np.quantile(y_pred_samples3, q=[0.250, 0.750], axis=1)

## 描画処理
# 描画領域の設定
plt.figure(figsize=(5, 4))
ax = plt.subplot()
# 観測値の散布図の描画
sns.scatterplot(ax=ax, data=data, x='X', y='Y', s=80, alpha=0.7, label='観測値')
# 予測値の中央値の折れ線グラフの描画
ax.plot(X_new, y_pred_median3, color='tab:red', label='予測値(中央値)')
# 予測値の95%CIの描画
ax.fill_between(X_new, y_pred_95ci3[0], y_pred_95ci3[1], color='tomato',
                 alpha=0.2, label='95%CI')
# 予測値の50%CIの描画
ax.fill_between(X_new, y_pred_50ci3[0], y_pred_50ci3[1], color='tomato',
                 alpha=0.5, label='50%CI')
# 修飾
ax.set(title='$y$ の観測値と予測値のプロット\nStudentのt分布',
       ylim=(-29, 79), yticks=range(-25, 76, 25), xticks=range(0, 12))
ax.legend(bbox_to_anchor=(1.41, 1))
ax.grid(lw=0.5);

【実行結果】

コーシー分布モデルとStudentのt分布モデルの比較

$${Y}$$の観測値と予測値のプロットを２つ横並びにして描画します。

## 描画横並び

# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 4), sharex=True, sharey=True)

# 観測値の散布図の描画
sns.scatterplot(ax=ax1, data=data, x='X', y='Y', s=80, alpha=0.7)
# 予測値の中央値の折れ線グラフの描画
ax1.plot(X_new, y_pred_median2, color='tab:red')
# 予測値の95%CIの描画
ax1.fill_between(X_new, y_pred_95ci2[0], y_pred_95ci2[1], color='tomato',
                 alpha=0.2)
# 予測値の50%CIの描画
ax1.fill_between(X_new, y_pred_50ci2[0], y_pred_50ci2[1], color='tomato',
                 alpha=0.5)
# 修飾
ax1.set(title='$y$ の観測値と予測値のプロット\nコーシー分布',
       ylim=(-29, 79), yticks=range(-25, 76, 25), xticks=range(0, 12))
ax1.grid(lw=0.5);


# 観測値の散布図の描画
sns.scatterplot(ax=ax2, data=data, x='X', y='Y', s=80, alpha=0.7, label='観測値')
# 予測値の中央値の折れ線グラフの描画
ax2.plot(X_new, y_pred_median3, color='tab:red', label='予測値(中央値)')
# 予測値の95%CIの描画
ax2.fill_between(X_new, y_pred_95ci3[0], y_pred_95ci3[1], color='tomato',
                 alpha=0.2, label='95%CI')
# 予測値の50%CIの描画
ax2.fill_between(X_new, y_pred_50ci3[0], y_pred_50ci3[1], color='tomato',
                 alpha=0.5, label='50%CI')
# 修飾
ax2.set(title='$y$ の観測値と予測値のプロット\nStudentのt分布',
       ylim=(-29, 79), yticks=range(-25, 76, 25), xticks=range(0, 12))
ax2.grid(lw=0.5)
ax2.legend(bbox_to_anchor=(1.5, 1))
plt.tight_layout();

【実行結果】
両方のモデルは、予測値（中央値、ベイズ信用区間）が$${X=3}$$の外れ値に引きずられていない様子が分かります。
コーシー分布のほうが95%CIの幅が狭い感じがします。

7.9 節は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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の教科書です。
よかったらぜひ、お試しくださいませ。

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

この記事が参加している募集

仕事について話そう

108,577件

この記事が気に入ったらサポートをしてみませんか？