見出し画像

StanとRでベイズ統計モデリングをPyMC Ver.5で写経～第７章「7.2 対数をとるか否か」

ネイピア DS

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

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

この記事は、テキスト第７章「回帰分析の悩みどころ」の7.2節「対数を取るか否か」の PyMC5写経 を取り扱います。
目的変数と説明変数の値を対数$${\log_{10}}$$で変換してモデリングします。

はじめに

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.2 対数をとるか否か

インポート

### インポート

# 数値・確率計算
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 matplotlib.cm as cm
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'

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

データの読み込み

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

### データの読み込み ◆data-rental.txt
# Y:2年間のトータル費用(万円), Area:物件の広さ(m²)

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

【実行結果】

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

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

【実行結果】

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

【実行結果】

テキスト図7.1の散布図を描画します。

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

# 描画領域の設定
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
# 左の散布図の描画（通常スケール）
sns.scatterplot(ax=ax1, data=data, x='Area', y='Y', s=80, alpha=0.3)
ax1.set(title='スケーリングなし')
ax1.grid(lw=0.5)
# 右の散布図の描画（両軸は対数スケール）
sns.scatterplot(ax=ax2, data=data, x='Area', y='Y', s=80, alpha=0.3)
ax2.set(xscale='log', xticks=(10, 20, 50, 100), xticklabels=(10, 20, 50, 100), 
        yscale='log', ylim=(100, 2500),
        yticks=(200, 500, 1000, 2000), yticklabels=(200, 500, 1000, 2000),
        title='両軸 対数スケール')
ax2.grid(lw=0.5)
plt.tight_layout();

【実行結果】

モデル式 7-1

データのスケールをそのまま使うモデルです。

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

# Yの予測分布に用いるAreaの値の設定
Area_new = np.linspace(10, 120, 50)

モデルの定義です。

### モデルの定義１ ◆モデル式7-1 通常スケール

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

    ### 事前分布
    b1 = pm.Uniform('b1', lower=-1000, upper=1000)
    b2 = pm.Uniform('b2', lower=-1000, upper=1000)
    sigma = pm.Uniform('sigma', lower=0, upper=1000)

    ### 線形予測子
    mu = pm.Deterministic('mu', b1 + b2 * Area, dims='data')
    
    ### 尤度関数
    obs = pm.Normal('obs', mu=mu, sigma=sigma, observed=Y, dims='data')

    ### 計算値
    yPred = pm.Normal('yPred', mu=mu, sigma=sigma, dims='data')
    yNew = pm.Normal('yNew', mu=b1 + b2 * AreaNew, 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.8,
                       nuts_sampler='numpyro', random_seed=1234)

【実行結果】省略

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

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

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

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

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

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

【実行結果】

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

### トレースプロットの表示
var_names = ['b1', 'b2', 'sigma', 'mu', 'yPred', '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 = ['b1', 'b2', 'sigma']
param_samples = idata1.posterior[vars].to_dataframe().reset_index(drop=True)
display(make_stats_df(param_samples).round(2))

【実行結果】

モデル式 7-2

目的変数と説明変数を対数変換するモデルです。

目的変数 Y、説明変数 Area、予測に用いる Area_new の値を対数$${\log_{10}}$$に変換します。

### 変数を対数スケールに変換
data['log_Y'] = np.log10(data['Y']).values
data['log_Area'] = np.log10(data['Area']).values
log_Area_new = np.log10(Area_new)

モデルの定義です。

### モデルの定義２ ◆モデル式7-2 対数スケール

with pm.Model() as model2:
    
    ### データ関連定義
    ## coordの定義
    model2.add_coord('data', values=data.index, mutable=True)
    ## dataの定義
    # 目的変数 Y
    Y = pm.ConstantData('Y', value=data['log_Y'].values, dims='data')
    # 説明変数 Area
    Area = pm.ConstantData('Area', value=data['log_Area'].values, dims='data')
    # 予測用の説明変数 AreaNew
    AreaNew = pm.ConstantData('AreaNew', value=log_Area_new, dims='dataNew')

    ### 事前分布
    b1 = pm.Uniform('b1', lower=-1000, upper=1000)
    b2 = pm.Uniform('b2', lower=-1000, upper=1000)
    sigma = pm.Uniform('sigma', lower=0, upper=1000)

    ### 線形予測子
    mu = pm.Deterministic('mu', b1 + b2 * Area, dims='data')
    
    ### 尤度関数
    obs = pm.Normal('obs', mu=mu, sigma=sigma, observed=Y, dims='data')

    ### 計算値
    yPred = pm.Normal('yPred', mu=mu, sigma=sigma, dims='data')
    yNew = pm.Normal('yNew', mu=b1 + b2 * AreaNew, sigma=sigma, dims='dataNew')

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

### モデルの表示
model2

【実行結果】

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

【実行結果】

MCMCを実行します。

### 事後分布からのサンプリング 20秒
with model2:
    idata2 = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.8,
                       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 = ['b1', 'b2', 'sigma', 'mu']
pm.summary(idata2, hdi_prob=0.95, var_names=var_names, round_to=3)

【実行結果】

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

### トレースプロットの表示
var_names = ['b1', 'b2', 'sigma', 'mu', 'yPred', 'yNew']
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 = ['b1', 'b2', 'sigma']
param_samples = idata2.posterior[vars].to_dataframe().reset_index(drop=True)
display(make_stats_df(param_samples).round(2))

【実行結果】

２つのモデルの比較

テキストの３種類のグラフを描画します。

■ 図7.2左 予測分布 モデル式7-1

### 目的変数の予測分布yNewの描画：通常スケール ◆図7.2左

## 描画用データの作成
# MCMCサンプリングデータからyNewを取り出し
y_pred_samples = idata1.posterior.yNew.stack(sample=('chain', 'draw')).data
# サンプリングデータの10%,50%,90%パーセンタイル点を算出してデータフレーム化
y_pred_df = pd.DataFrame(
    np.quantile(y_pred_samples, q=[0.1, 0.25, 0.5, 0.75, 0.9], axis=1).T,
    columns=['10%', '25%', 'median', '75%', '90%'])

## 描画処理
# 観測値の散布図の描画
sns.scatterplot(data=data, x='Area', y='Y', s=80, alpha=0.5, label='観測値')
# yの予測分布の中央値を描画
plt.plot(Area_new, y_pred_df['median'], color='tomato', label='中央値')
# yの予測分布の80%CIを描画
plt.fill_between(Area_new, y_pred_df['10%'].values, y_pred_df['90%'].values,
                 color='tomato', alpha=0.2, label='80%CI')
# yの予測分布の50%CIを描画
plt.fill_between(Area_new, y_pred_df['25%'].values, y_pred_df['75%'].values,
                 color='tomato', alpha=0.3, label='50%CI')
# 修飾
plt.title('$Y$のベイズ予測分布(通常スケール)')
plt.grid(lw=0.5)
plt.legend();

【実行結果】

■ 図7.2右 予測分布 モデル式7-2

### 目的変数の予測分布yNewの描画：対数スケール ◆図7.2右

## 描画用データの作成
# MCMCサンプリングデータからyNewを取り出し（10のy_pred乗）
y_pred_samples_log = np.power(10, idata2.posterior.yNew.
                              stack(sample=('chain', 'draw')).data)
# サンプリングデータの10%,50%,90%パーセンタイル点を算出してデータフレーム化
y_pred_df_log = pd.DataFrame(
    np.quantile(y_pred_samples_log, q=[0.1, 0.25, 0.5, 0.75, 0.9], axis=1).T,
    columns=['10%', '25%', 'median', '75%', '90%'])

## 描画処理
# 観測値の散布図の描画
sns.scatterplot(data=data, x='Area', y='Y', s=80, alpha=0.5, label='観測値')
# yの予測分布の中央値を描画
plt.plot(Area_new, y_pred_df_log['median'], color='tomato', label='中央値')
# yの予測分布の80%CIを描画
plt.fill_between(Area_new,
                 y_pred_df_log['10%'].values, y_pred_df_log['90%'].values,
                 color='tomato', alpha=0.2, label='80%CI')
# yの予測分布の50%CIを描画
plt.fill_between(Area_new,
                 y_pred_df_log['25%'].values, y_pred_df_log['75%'].values,
                 color='tomato', alpha=0.3, label='50%CI')
# 修飾
plt.title('$Y$のベイズ予測分布(対数 $\log_{10}$)')
plt.grid(lw=0.5)
plt.legend();

【実行結果】

■ 図7.3左 観測値と予測値のプロット モデル式7-1

### 観測値と予測値のプロット：通常スケール ◆図7.3左

## 描画用データの作成 yPredの中央値と80%区間を算出
# MCMCサンプリングデータからyPredを取り出し
y_pred_samples2 = idata1.posterior.yPred.stack(sample=('chain', 'draw')).data
# サンプリングデータの10%,50%,90%パーセンタイル点を算出してデータフレーム化
y_pred_df2 = pd.DataFrame(
    np.quantile(y_pred_samples2, q=[0.1, 0.5, 0.9], axis=1).T,
    columns=['10%', 'median', '90%'])
y_pred_df2 = pd.concat([data, y_pred_df2], axis=1)
# 中央値と10%点の差、90%点と中央値の差を算出: errorbarで利用
y_pred_df2['err_lower'] = y_pred_df2['median'] - y_pred_df2['10%'] 
y_pred_df2['err_upper'] = y_pred_df2['90%'] - y_pred_df2['median']

## 描画処理
# 描画領域の指定
plt.figure(figsize=(6, 6))
ax = plt.subplot()
# 描画（エラーバー付き散布図）
ax.errorbar(y_pred_df2['Y'], y_pred_df2['median'],
            yerr=[y_pred_df2['err_lower'], y_pred_df2['err_upper']],
            color='tab:blue', alpha=0.5, marker='o', ms=10, linestyle='none')
# 赤い対角線の描画
ax.plot([-100, 1800], [-100, 1800], color='red', ls='--')
# 修飾
ax.set(xlabel='Observed: $Y$の観測値', ylabel='Predicted: $Y$の予測値(中央値)',
       title='$Y$ の観測値と予測値（中央値）のプロット\n80%区間, 通常スケール')
ax.grid(lw=0.5);

【実行結果】

■ 図7.3右 観測値と予測値のプロット モデル式7-2

### 観測値と予測値のプロット：対数スケール ◆図7.3右

## 描画用データの作成 yPredの中央値と80%区間を算出
# MCMCサンプリングデータからyPredを取り出し（10^y_pred乗）
y_pred_samples2_log = np.power(10, idata2.posterior.yPred.
                                   stack(sample=('chain', 'draw')).data)
# サンプリングデータの10%,50%,90%パーセンタイル点を算出してデータフレーム化
y_pred_df2_log = pd.DataFrame(
    np.quantile(y_pred_samples2_log, q=[0.1, 0.5, 0.9], axis=1).T,
    columns=['10%', 'median', '90%'])
y_pred_df2_log = pd.concat([data, y_pred_df2_log], axis=1)
# 中央値と10%点の差、90%点と中央値の差を算出: errorbarで利用
y_pred_df2_log['err_lower'] = y_pred_df2_log['median'] - y_pred_df2_log['10%'] 
y_pred_df2_log['err_upper'] = y_pred_df2_log['90%'] - y_pred_df2_log['median']

## 描画処理
# 描画領域の指定
plt.figure(figsize=(6, 6))
ax = plt.subplot()
# 描画（エラーバー付き散布図）
ax.errorbar(y_pred_df2_log['Y'], y_pred_df2_log['median'],
            yerr=[y_pred_df2_log['err_lower'], y_pred_df2_log['err_upper']],
            color='tab:blue', alpha=0.5, marker='o', ms=10, linestyle='none')
# 赤い対角線の描画
ax.plot([-100, 1800], [-100, 1800], color='red', ls='--')
# 修飾
ax.set(xlabel='Observed: $Y$の観測値', ylabel='Predicted: $Y$の予測値(中央値)',
       title='$Y$ の観測値と予測値（中央値）のプロット\n80%区間, '
             '対数 $\log_{10}$スケール')
ax.grid(lw=0.5);

【実行結果】

■ 図7.4左 推定されたノイズの分布 モデル式7-1

### 推定されたノイズの分布の描画 ◆図7.4左

## 描画用データの作成 yPredの中央値と80%区間を算出
# MCMCサンプリングデータからyPredを取り出し
mu_samples = idata1.posterior.mu.stack(sample=('chain', 'draw')).data
# εの算出： 観測値Y - muの事後分布サンプリングデータ
epsilon_samples = np.array([data.loc[i, 'Y'] - mu_sample
                            for i, mu_sample in enumerate(mu_samples)])

## εのMAP推定値の算出
# εのMAP推定値リストの初期化
map_plot = []
# 物件ごとのε（ε[n]）のKDEの最大値の取得
for i, epsilon in enumerate(epsilon_samples):
    # x軸の値を設定
    x_ticks = np.linspace(epsilon.min(), epsilon.max(), 1001)
    # KDEプロットの描画
    kde_plot = sns.kdeplot(epsilon, color=None)
    # KDEプロットで描画した線から、x軸,y軸の値を取得
    kde_data = kde_plot.get_lines()[i].get_data()
    # KDEプロットのyの最大値のインデックスを取得
    kde_max_idx = kde_data[1].argmax()
    # MAP推定値の格納
    map_plot.append(kde_data[0][kde_max_idx])
plt.close()
map_plot = np.array(map_plot)

## σのMAP推定値の算出
sigma_samples = idata1.posterior.sigma.stack(sample=('chain', 'draw')).data
sigma_plot = sns.kdeplot(sigma_samples)
sigma_kde_data = sigma_plot.get_lines()[0].get_data()
sigma_map = sigma_kde_data[0][sigma_kde_data[1].argmax()]
plt.close()

## 描画用データの作成
# 正規分布のx軸の値を設定
xvals = np.linspace(-400, 600, 1001)
# 正規分布の確率密度関数の算出
yvals = stats.norm.pdf(xvals, loc=0, scale=sigma_map)

## 描画処理 twinxを用いて、度数と密度の2つのyに対応
# 描画領域の設定
fig, ax1 = plt.subplots()
ax2 = ax1.twinx()            # y軸を左と右の２つにする
# ax1: εのMAP推定値のヒストグラムの描画
sns.histplot(map_plot, bins=20, color='tab:blue', ec='white', alpha=0.7,
             label='$\epsilon$のMAP推定値',
             ax=ax1)
# ax2: εのMAP推定値のKDEプロットの描画
sns.kdeplot(map_plot, color='tab:blue',
            label='$\epsilon$のMAP推定値のKDE',
            ax=ax2)
# ax2: σのMAP推定値を用いた正規分布の確率密度関数の描画
ax2.plot(xvals, yvals, color='tab:orange', label='正規分布(0, $\sigma$)')
# ax1とax2の凡例を統合
h1, l1 = ax1.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
ax1.legend(h1+h2, l1+l2, bbox_to_anchor=(1.55, 1))
# 修飾
ax1.set(xlabel='value')
ax1.set_title('$\epsilon[n]$ のMAP推定値のヒストグラムとKDE,\n'
              '$\sigma$ のMAP推定値を用いた正規分布の確率密度関数\n'
              '通常スケール')
ax1.grid(lw=0.5);

【実行結果】

■ 図7.4右 推定されたノイズの分布 モデル式7-2

### 推定されたノイズの分布の描画 ◆図7.4右

## 描画用データの作成 yPredの中央値と80%区間を算出
# MCMCサンプリングデータからyPredを取り出し
mu_samples = idata2.posterior.mu.stack(sample=('chain', 'draw')).data
# εの算出： log10(観測値Y) - muの事後分布サンプリングデータ
epsilon_samples = np.array([np.log10(data.loc[i, 'Y']) - mu_sample
                            for i, mu_sample in enumerate(mu_samples)])

## εのMAP推定値の算出
# εのMAP推定値リストの初期化
map_plot = []
# 物件ごとのε（ε[n]）のKDEの最大値の取得
for i, epsilon in enumerate(epsilon_samples):
    # x軸の値を設定
    x_ticks = np.linspace(epsilon.min(), epsilon.max(), 1001)
    # KDEプロットの描画
    kde_plot = sns.kdeplot(epsilon, color=None)
    # KDEプロットで描画した線から、x軸,y軸の値を取得
    kde_data = kde_plot.get_lines()[i].get_data()
    # KDEプロットのyの最大値のインデックスを取得
    kde_max_idx = kde_data[1].argmax()
    # MAP推定値の格納
    map_plot.append(kde_data[0][kde_max_idx])
plt.close()
map_plot = np.array(map_plot)

## σのMAP推定値の算出
sigma_samples = idata2.posterior.sigma.stack(sample=('chain', 'draw')).data
sigma_plot = sns.kdeplot(sigma_samples)
sigma_kde_data = sigma_plot.get_lines()[0].get_data()
sigma_map = sigma_kde_data[0][sigma_kde_data[1].argmax()]
plt.close()

## 描画用データの作成
# 正規分布のx軸の値を設定
xvals = np.linspace(-0.5, 0.5, 1001)
# 正規分布の確率密度関数の算出
yvals = stats.norm.pdf(xvals, loc=0, scale=sigma_map)

## 描画処理 twinxを用いて、度数と密度の2つのyに対応
# 描画領域の設定
fig, ax1 = plt.subplots()
ax2 = ax1.twinx()            # y軸を左と右の２つにする
ax3 = ax1.twinx()            # y軸を左と右の２つにする
# ax1: εのMAP推定値のヒストグラムの描画
sns.histplot(map_plot, bins=20, color='tab:blue', ec='white', alpha=0.7,
             label='$\epsilon$のMAP推定値',
             ax=ax1)
# ax2: εのMAP推定値のKDEプロットの描画
sns.kdeplot(map_plot, color='tab:blue',
            label='$\epsilon$のMAP推定値のKDE',
            ax=ax2)
ax2.axis('off')
# ax3: σのMAP推定値を用いた正規分布の確率密度関数の描画
ax3.plot(xvals, yvals, color='tab:orange', label='正規分布(0, $\sigma$)')
ax3.axis('off')
# ax1～ax3の凡例を統合
h1, l1 = ax1.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
h3, l3 = ax3.get_legend_handles_labels()
ax1.legend(h1+h2+h3, l1+l2+l3, bbox_to_anchor=(1.4, 1))
# 修飾
ax1.set(xlabel='value')
ax1.set_title('$\epsilon[n]$ のMAP推定値のヒストグラムとKDE,\n'
              '$\sigma$ のMAP推定値を用いた正規分布の確率密度関数\n'
              '対数 $\log_{10}$スケール')
ax1.grid(lw=0.5);

【実行結果】

7.2 節は以上です。

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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

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

いいなと思ったら応援しよう！

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

#仕事について話そう

115,662件

#うちの積読を紹介する

1,892件