「ベイズ統計モデリングによるデータ分析入門」をPythonで写経 ～ Vol.2 データの要約

書籍の著者　馬場真哉 先生

この記事は、書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」第２部第２章「データの要約」の Python 写経活動記録です。

この章の記述統計手法は、MCMCサンプルからパラメータを推定する際によく用いられるものです。
各種手法を Python で動かして体感し、見慣れておきましょう！

では書籍を開いてベイズ統計モデリングの旅に出かけましょう🚀

はじめに

このブログシリーズは書籍「RとStanではじめるベイズ統計モデリングによるデータ分析入門」（講談社、「テキスト」と呼びます）の Python 写経です。

テキストの紹介と引用表記はリンク先の記事に掲載しています。

準備

■ 記事の範囲
この記事はテキスト第２部第２章の以下の節を取り扱います。

■ 利用データ
テキスト・サポートサイトのデータファイルを引用しています。
Jupyter Notebook ファイルと同一フォルダ内の「data」フォルダにデータファイルを格納しています。

■ ライブラリのインポート
Jupyter Notebook 形式でコードを記述します。
この記事で用いるライブラリをインポートします。

# インポート

# 数値計算
import numpy as np
import pandas as pd

# 統計
import scipy.stats as stats
import statsmodels.api as sm                 # コレログラム
import statsmodels.tsa.api as tsa            # 自己共分散、自己相関係数

# 機械学習
from sklearn.neighbors import KernelDensity  # カーネル密度推定

# 可視化
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'  # または import japanize_matplotlib 等

第２章 データの要約

度数・度数分布・ヒストグラム

データの csv ファイルを pandas のデータフレーム形式で変数 fish に読み込みます。
read_csv 関数を利用します。

# p.93 データの読み込み
fish = pd.read_csv('./data/2-2-1-fish.csv')
print('fish.shape: ', fish.shape)
fish.head(3)

【実行結果】
変数 length（体長）に 400 個のデータが含まれています。

ヒストグラムを描画します。
pandas データフレームの plot メソッドを利用します。
図 2.2.1 に相当します。

# p.94 ヒストグラム ※pandasのplotメソッド利用
fish.plot.hist(edgecolor='white', alpha=0.7, xlabel='体長 [cm]', ylabel='頻度');

【実行結果】
ベル型の正規分布に似た形状をしています。

🔷🔷🔷

カーネル密度推定

ベイズ統計モデリングの文脈でカーネル密度推定は、MCMCサンプルを用いたパラメータの事後分布等の可視化（plot_trace や plot_density）などで利用されます。

fish データに関するカーネル密度推定を実行して、図 2.2.5 相当の曲線を描きます。
scipy.stats の gaussian_kde クラスを利用します。
インスタンス作成 ⇒ 推定 の２ステップです。

# p.96 図2.2.5 カーネル密度推定 ※scipy.statsのガウスカーネル密度推定を利用

## 設定
# x軸の値の設定
x_val = np.linspace(2, 17, 1001)

## カーネル密度推定
# インスタンス作成
kde = stats.gaussian_kde(fish['length'])
# x_valに関するカーネル密度推定
kde_val = kde(x_val)

## 描画
# カーネル密度推定した曲線の描画
plt.figure(figsize=(8, 3))
plt.plot(x_val, kde_val)
# 修飾
plt.title('カーネル密度推定')
plt.xlabel('体長', fontsize=12)
plt.ylabel('Density', fontsize=12)
plt.show()

【実行結果】
先程のヒストグラムがやや滑らかな曲線に変換されました。

カーネル密度推定の形状に影響を与える引数が「バンド幅」です。
「３つのバンド幅調整値」を比較する図 2.2.6 相当の描画に取り組みます。
３つのライブラリを使ってみて、比べてみましょう。

1️⃣ scipy.stats の gaussian_kde クラス
先ほどと同じクラスです。
for 文でバンド幅の調整値を切り替えて、カーネル密度推定の算出と曲線の描画を行います。

# p.97 図2.2.6 バンド幅の変更 ※scipy利用

## 設定
# バンド
band_adjusts = [1, 0.25, 4]
# 線の色
colors = ['black', 'tab:orange', 'tab:blue']

## カーネル密度推定の実行
# インスタンス作成
kde_sp = stats.gaussian_kde(fish['length'])
# 標準のバンド幅の取得
band_std_sp = kde_sp.factor

## 描画
# 描画領域の設定
plt.figure(figsize=(8, 3))
# ３つのバンド幅調整値ごとにKDEプロット描画を繰り返し処理
for adjust, color in zip(band_adjusts, colors):
    # バンド幅をadjust倍して設定
    kde_sp.set_bandwidth(band_std_sp * adjust)
    # x_valに関するカーネル密度推定
    kde_sp_val = kde_sp(x_val)
    # KDEプロットの描画
    plt.plot(x_val, kde_sp_val, label=f'バンド幅{adjust}倍', color=color)
# 修飾
plt.title('カーネル密度推定 scipy')
plt.xlabel('体長', fontsize=12)
plt.ylabel('Density', fontsize=12)
plt.legend()
plt.show()

【実行結果】
バンド幅によって曲線の形状が全然違ってきますね。

2️⃣ seaborn の kde_plot
seaborn は 可視化ライブラリです。
kde_plot で自動的にカーネル密度推定と可視化を同時実行してくれます。
裏側では scipy の gaussian_kde クラスが動いているそうです。

# p.97 図2.2.6 バンド幅の変更 ※seaborn利用（バックエンドでscipy利用）

## 描画領域の設定
plt.figure(figsize=(8, 3))
# ３つのバンド幅調整値ごとにKDEプロット描画を繰り返し処理
for adjust, color in zip(band_adjusts, colors):
    # KDEプロットの描画
    sns.kdeplot(
        data=fish['length'],        # データ
        bw_adjust=adjust,           # バンド幅調整値
        color=color,                # 線の色
        label=f'バンド幅{adjust}倍'  # 凡例のラベル 
    )
# 修飾
plt.title('カーネル密度推定 seaborn')
plt.xlabel('体長', fontsize=12)
plt.ylabel('Density', fontsize=12)
plt.xlim(2, 17)
plt.legend()
plt.show()

【実行結果】
可視化だけを行いたい場合には、seaborn を使うのが良さそうですね。

3️⃣ scikit-learn の KernelDensity クラス
scikit-learn は機械学習用のライブラリです。

# p.97 図2.2.6 バンド幅の変更 ※scikit-learn利用

## カーネル密度推定の実行
# カーネル密度推定インスタンスの作成と学習
kde_sk = KernelDensity(kernel='gaussian', bandwidth='scott').fit(fish.values)
# 標準のバンド幅の取得 ※scipyのバンド幅計算と一致するようにデータの標準偏差を乗じる
band_std_sk = kde_sk.bandwidth_ * fish['length'].std(ddof=1)

## 描画
# 描画領域の設定
plt.figure(figsize=(8, 3))
# ３つのバンド幅調整値ごとにKDEプロット描画を繰り返し処理
for adjust, color in zip(band_adjusts, colors):
    # 調整後のバンド幅でカーネル密度推定の実行
    kde_sk_tmp = KernelDensity(
        kernel='gaussian',            # カーネルの種類
        bandwidth=band_std_sk *adjust  # バンド幅
    ).fit(fish.values)
    # カーネル密度推定値の算出（対数スケール）
    log_density = kde_sk_tmp.score_samples(x_val.reshape(-1, 1))
    # KDEプロットの描画
    plt.plot(x_val, np.exp(log_density), label=f'バンド幅{adjust}倍', color=color)
# 修飾
plt.title('カーネル密度推定 scikit-learn')
plt.xlabel('体長', fontsize=12)
plt.ylabel('Density', fontsize=12)
plt.legend()
plt.show()

【実行結果】

【コードの補足】
scipy のカーネル密度推定結果と合わせる目的で次のようにしています。

• 標準のバンド幅の算出
  以下を行って標準のバンド幅を計算しています。

  • KernelDensity の引数 bandwidth でバンド幅推定法「scott」を指定

  • 最初の KernelDensity で推定したバンド幅にデータの標準偏差を乗算
    （band_std_sk = kde_sk.bandwidth_ * fish['length'].std(ddof=1)）

• 最終的なバンド幅の算出
  描画時に「バンド幅＝標準のバンド幅 $${\times}$$ バンド幅調整値」を用いて再度、カーネル密度推定を行っています。

🔷🔷🔷

算術平均

いわゆる平均です。
pandas データフレーム の mean メソッドを利用します。
データフレーム変数 fish の後ろに mean() を指定します。

# p.97 算術平均
fish.mean()

【実行結果】

🔷🔷🔷

中央値・四分位点・パーセント点

0 から 1000 の等差級列（連続する整数値）を作成します。
numpy の arange() 関数を用いて、結果は numpy 配列で出力されます。

# p.98 等差数列の作成
suuretu = np.arange(1001)
suuretu

【実行結果】

この配列の要素数（長さ）を算出します。
python 組み込み関数の len を使います。

# p.98 数列の長さ
len(suuretu)

【実行結果】
1001 個の要素があります。

この配列の中央値を算出します。
numpy の median() 関数を利用します。

# p.98 数列の中央値
np.median(suuretu)

【実行結果】
ちょうど真ん中の数値は 500 です。

中央値は 50% 点です。
numpy の quantile （分位数）関数で配列の 50% 点（q=0.5）を求めます。

# p.98 数列の中央値 quantile
np.quantile(suuretu, q=0.5)

【実行結果】

続いて第１四分位数（25% 点）と第３四分位数（75% 点）を求めます。
numpy の quantile 関数の引数 q= [0.25, 0.75] を与えます。

# p.99 数列の25%点、75%点
np.quantile(suuretu, q=[0.25, 0.75])

【実行結果】

続いて 95% 区間（2.5% 点と 97.5% 点）を求めます。

# p.99 数列の2.5%点、97.5%点
np.quantile(suuretu, q=[0.025, 0.975])

【実行結果】

ベイズ統計モデリングでは、パラメータの事後分布から生成した「MCMC サンプル」に関して、平均値：事後平均、中央値：事後中央値、 95% 区間：95% ベイズ信用区間といった統計量をよく見ます！

🔷🔷🔷

共分散とピアソンの積率相関係数

新しいデータの csv ファイルを pandas の データフレーム birds に読み込みます。

# p.99 ファイルの読み込み

birds = pd.read_csv('./data/2-1-1-birds.csv')
print('birds.shape: ', birds.shape)
birds.head(3)

【実行結果】
標本サイズは８、列（変数）は３つです。

体の大きさ body_length と 羽の長さ feather_length の相関係数を求めます。
pandas データフレームの corr メソッドを利用します。
引数 numeric_only=True とすることで、数値型の変数だけを与えられます。

# p.99 体の大きさと羽の大きさの相関係数
birds.corr(numeric_only=True)

【実行結果】
相関係数は 0.996 。ものすごく強い正の相関がありますね…

🔷🔷🔷

自己共分散・自己相関係数・コレログラム

新しいデータの csv ファイルを pandas の データフレーム nile に読み込みます。
このファイルは WEB サイトから直接取得できます。

# p.100 ナイル川流量データの読み込み WebサイトのCSVファイルを直接読み込む

# ナイル川流量データのURL
url = 'https://vincentarelbundock.github.io/Rdatasets/csv/datasets/Nile.csv'
# CSVファイルの読み込み
nile = pd.read_csv(url, usecols=[1, 2], index_col=0)
# 結果の表示
print('nile.shape: ', nile.shape)
nile

【実行結果】
1871～1970年のナイル川の流量 value です。

標本自己共分散を算出します。
python の 統計ライブラリ statsmodels の acovf 関数を利用します。
テキストの max.lag = 5 を 引数 nlag で与えます。

# p.100 標本自己共分散 statsmodels利用

# ラグの最大値
nlag = 5
# 標本自己共分散の算出
nile_acov = tsa.acovf(nile, nlag=nlag)
# 結果をデータフレームで表示
pd.DataFrame(nile_acov, index=range(nlag+1), columns=['自己共分散']).T.round(1)

【実行結果】
ラグ０～５の自己共分散です。

ラグは期間のズレを指します。
ラグ１は１期ずれ、ラグ０はズレなしです。
「自己」共分散は、同一データ、つまり「自分自身」の過去期間との共分散です。
ラグ１の場合は、１期ずらした自己データとの共分散です。

続いて標本自己相関係数を算出します。
statsmodels の acf 関数を利用します。

# p.101 標本自己相関係数 statsmodels利用

# ラグの最大値
nlag = 5
# 標本自己相関係数の算出
nile_acf = tsa.acf(nile, nlags=nlag)
# 結果をデータフレームで表示
pd.DataFrame(nile_acf, index=range(nlag+1), columns=['自己共分散']).T.round(4)

【実行結果】

自己相関係数 ACF を描画しましょう。
「コレログラム」と呼ばれるチャートです。
statsmodels の plot_acf 関数を利用します。
引数 lag = 20（ラグ20まで計算）などを与えています。

# p.101 コレログラムの描画 statsmodels利用
sm.graphics.tsa.plot_acf(nile, lags=20, auto_ylims=True);

【実行結果】
青塗り領域（95% 信頼区間）を飛び抜けているラグ３までが有意な感じです。

ベイズ統計モデリングでは、MCMC サンプルに自己相関が生じていないかどうかを確認することがあります！

今回の記事は以上です。
楽しかったですね！

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

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

１．のんびり統計

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

２．統計・データ分析とつながる

シリーズ「統計・データ分析とつながる」は、統計・データ分析との「つながり」を発掘して、コラム風に仕立てたブログシリーズです。
生成 AI の力を借りながら、統計・データ分析の入り口をイメージして、自由気ままに書きました。
たとえば…
・日常生活と統計のつながり
・統計検定２級からその先へのつながり
気楽にお読みいただけたら嬉しいです🍀

３．実験！たのしいベイズモデリング１＆２を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の教科書です。
よかったらぜひ、お試しくださいませ。

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