「数値シミュレーションで読み解く統計のしくみ」をPythonで写経～第１章「本書のねらい」

ネイピア DS

2024年6月15日 11:46

第１章「本書のねらい」

書籍の著者　小杉考司先生、紀ノ定保礼先生、清水裕士先生

この記事は、書籍「数値シミュレーションで読み解く統計のしくみ」第１章「本書のねらい」の Python写経活動 を取り扱います。

今回は「シミュレーションのデモンストレーション」を写経します！
では書籍を開いてシミュレーションの旅に出発です🚀

はじめに

書籍「数値シミュレーションで読み解く統計のしくみ」のご紹介

このシリーズは書籍「数値シミュレーションで読み解く統計のしくみ〜Rでためしてわかる心理統計」（技術評論社、「テキスト」と呼びます）の Python写経です。

テキストは、2023年9月に発売され、副題「Rでためしてわかる心理統計」のとおり、統計処理に定評のある R の具体的なコードを用いて、心理統計の理解に役立つ数値シミュレーションを実践する素晴らしい書籍です。

大まかな構成は以下のとおりです。

R の練習
シミュレーションのベースになる「乱数生成」の練習
母数推定のシミュレーション
統計的検定のエラー確率のシミュレーション
サンプルサイズ設計のシミュレーション
回帰分析のシミュレーション

とにかくサンプルコードの量が多いです！
確率・統計分野の主だったライブラリを余す所なく満喫できるのです！
（個人の見解です🍀）
コードをタンマリと動かせてニンマリなのです。
（個人の見解です🍀）

Python写経って？

R のサンプルコードを Python 化 します！
統計を学ぶ方の中には Python でやってみたいという人も居ることでしょう。
お薦めしたいこのテキストを Python で楽しむ仲間が増えたなら嬉しいです✨️

R は Python に比べて統計処理がとても充実しているので、もしかすると Python で実現できないことが現れるかもしれません。
それはそれで Python の限界（というよりも私の翻訳能力の限界）が明らかになって、楽しめそうです！

引用表記

この記事は、出典に記載の書籍に掲載された文章及びコードを引用し、適宜、掲載文章とコードを改変して書いています。
【出典】
「数値シミュレーションで読み解く統計のしくみ〜Rでためしてわかる心理統計」初版第１刷、著者小杉考司・紀ノ定保礼・清水裕士、技術評論社

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

第１章本書のねらい

テキストは、「確率変数」を扱う統計分析が、母集団や標本といった具体的な数字を超えた「抽象的な実体」を扱っており、この抽象的な実体をイメージできないことで解析結果の解釈を誤るリスクを警鐘しています。
この抽象さに具体性を与える方法としてシミュレーションを推奨しています。
第１章ではRを用いて「シミュレーションの例を紹介」しています。
R の文法はひとまず置いておき、デモンストレーションです。

もちろんこの記事は R ではなく Python で描きます。
Jupyter Notebook 形式（拡張子 .ipynb）でPythonコードを書きます。

この記事で用いるライブラリをインポートします。

### インポート
# 数値計算
import numpy as np
# 確率・統計
import scipy.stats as stats
# 描画
import matplotlib.pyplot as plt

シミュレーション例１

母集団の相関係数が 0.5 のときに、データから得られる標本相関係数が 0.5 からどの程度変動するかシミュレーションして可視化します。

まずはサンプルサイズ$${n=25}$$のケースです。
２つの乱数を25個ずつ作成して相関係数を算出する処理を1000回繰り返します。
標本相関係数を1000個取得しています。
テキストの図1.2を描画します。

### 4ページ 母相関係数0.5によるサンプルサイズn=25の乱数生成 5ページ図1.2

## 設定と準備
rho = 0.5     # 母相関係数
n = 25        # サンプルサイズ
iter = 1000   # 反復回数

# 結果を格納するオブジェクト
r = np.zeros(iter)

## シミュレーション
# 乱数生成器（乱数シード）の設定
rng = np.random.default_rng(seed=123)
for i in range(iter):
    # Y1の作成 標準正規分布乱数
    Y1 = rng.normal(size=n, loc=0, scale=1)
    # Y2の作成 Y1に母相関係数rhoを乗じて正規分布のノイズを加算
    Y2 = Y1 * rho + rng.normal(size=n, loc=0, scale=(1 - rho**2)**0.5)
    # Y1とY2の標本相関係数を算出
    r[i] = np.corrcoef(Y1, Y2)[0, 1]

## 結果
plt.hist(r, edgecolor='white')
plt.xlim(-0.2, 1)
plt.title('Histogram of r')
plt.xlabel('r')
plt.ylabel('Frequency');

【実行結果】
テキストの解説通り、母相関係数 0.5 を設定したにもかかわらず、シミュレーションでは負の相関係数が発生しています。

【コード補足】
１．乱数生成
numpy の Random Generator を利用しています。
この写経シリーズでは乱数生成になるべくnumpy の Random Generator を用いるようにします。
rng = np.random.default_rng(seed=123) で、乱数シードに123を設定した numpy 乱数生成器を作成します。
以後の乱数生成時には、rng.xxx のように、乱数生成器 rng から乱数生成を呼び出します。
たとえば、rng.normal(loc=平均, scale=標準偏差, size=乱数生成数) で正規分布乱数を生成します。

２．相関係数の算出
numpy の corrcoef() を利用しています。
np.corrcoef(変数1, 変数2) で変数1 と変数2 の相関係数を算出できます。

続いてサンプルサイズ$${n=100}$$のケースです。
テキストの図1.3を描画します。

### 5ページ 母相関係数0.5によるサンプルサイズn=100の乱数生成 図1.3

## 設定と準備
rho = 0.5     # 母相関の設定
n = 100       # サンプルサイズ
iter = 1000   # 反復回数

# 結果を格納するオブジェクト
r = np.zeros(iter)

## シミュレーション
rng = np.random.default_rng(seed=123)
for i in range(iter):
    Y1 = rng.normal(size=n, loc=0, scale=1)
    Y2 = Y1 * rho + rng.normal(size=n, loc=0, scale=(1 - rho**2)**0.5)
    r[i] = np.corrcoef(Y1, Y2)[0, 1]

## 結果
plt.hist(r, edgecolor='white')
plt.xlim(-0.2, 1)
plt.title('Histogram of r')
plt.xlabel('r')
plt.ylabel('Frequency');

【実行結果】
負の相関係数は無くなりましたが、まだまだ 0.2 ～ 0.7 の幅があります。
標本から得られる統計量の解釈って本来難しいものなのかもです・・・

シミュレーション例２

テキストは無相関検定のシミュレーションに移ります。
無相関とは相関係数が０のことです。
無相関検定の帰無仮説は「２つのデータは無相関である」です。

無相関の乱数データを生成して無相関検定を行い、p値を取得します。
無相関になるように生成したデータに対する無相関検定で、有意水準$${\alpha=0.05}$$を下回る「有意になる確率」のシミュレーションです。
誤って帰無仮説を棄却する確率をシミュレーションで割り出すのだと思います。

まずはサンプルサイズ$${n=25}$$のケースです。

### 6ページ 無相関検定 サンプルサイズ25, 10000回検定

## 設定と準備
rho = 0       # 母相関の設定
n = 25        # サンプルサイズ
iter = 10000  # 反復回数

# 結果を格納するオブジェクト
p = np.zeros(iter)

## シミュレーション
rng = np.random.default_rng(seed=123)         # 乱数生成器（乱数シード）
for i in range(iter):
    Y1 = rng.normal(size=n, loc=0, scale=1)
    Y2 = Y1 * rho + rng.normal(size=n, loc=0, scale=(1 - rho**2)**0.5)
    p[i] = stats.pearsonr(Y1, Y2).pvalue      # 無相関検定のp値、両側検定

## 結果
(p < 0.05).mean()

【実行結果】
有意水準$${\alpha=0.05}$$を下回った割合（確率）は$${0.493}$$。
有意水準と同レベルです！

【コード補足】無相関検定
scipy.stats の pearsonr() を利用しています。
scipy.stats を stats の名前でインポートしていますので、stats.pearsonr() で呼び出します。
stats.pearsonr(変数1, 変数2) で変数1 と変数2 の無相関検定を行えます。
stats.pearsonr(変数1, 変数2).pvalue で無相関検定のp値を取得できます。

続いてサンプルサイズ$${n=2500}$$のケースです。

### 6ページ 無相関検定 サンプルサイズ2500, 10000回検定

## 設定と準備
rho = 0       # 母相関の設定
n = 2500      # サンプルサイズ
iter = 10000  # 反復回数

# 結果を格納するオブジェクト
p = np.zeros(iter)

## シミュレーション
rng = np.random.default_rng(seed=123)
for i in range(iter):
    Y1 = rng.normal(size=n, loc=0, scale=1)
    Y2 = Y1 * rho + rng.normal(size=n, loc=0, scale=(1 - rho**2)**0.5)
    p[i] = stats.pearsonr(Y1, Y2).pvalue

## 結果
(p < 0.05).mean()

【実行結果】
5%近くに収まり、有意水準と同レベルであることが分かりました。

有意水準の意義になるほどです！
テキストの言葉をお借りします。

このように、適切な検定が行われれば、サンプルサイズによらず有意になる確率は設定した有意水準に一致するようになります。

テキストより引用

シミュレーション例３

テキストは誤った統計分析に関するシミュレーションを実行して、問題点を明瞭に可視化しています。
「もう少しで有意水準をクリアできそうだから、ちょっとデータを増やしてみよう」というｎ水増し問題のシミュレーションです。
４～５分ほど処理時間がかかりました。

### 7ページ n増し問題 無相関検定 サンプルサイズ25, 10000回検定 4分30秒

## 設定と準備
rho = 0       # 母相関の設定
n = 25        # サンプルサイズ
iter = 10000  # 反復回数
alpha = 0.05  # 有意水準

# 結果を格納するオブジェクト
p = np.zeros(iter)

## シミュレーション
rng = np.random.default_rng(seed=123)           # 乱数生成器（乱数シード）
for i in range(iter):
    # 最初のデータ
    Y1 = rng.normal(size=n, loc=0, scale=1)
    Y2 = Y1*rho + rng.normal(size=n, loc=0, scale=(1 - rho**2)**0.5)
    p[i] = stats.pearsonr(Y1, Y2).pvalue        # 無相関検定のp値、両側検定
    # データ追加
    count = 0
    ## p値が5%を下回るか、データが当初の倍になるまで増やし続ける
    while p[i] >= alpha and count < n * 2:
        # 有意ではなかったとき、それぞれの変数に1つデータを追加
        Y1_add = rng.normal(size=1, loc=0, scale=1)
        Y1 = np.append(Y1, Y1_add)
        Y2_add = Y1_add*rho + rng.normal(size=1, loc=0, scale=(1 - rho**2)**0.5)
        Y2 = np.append(Y2, Y2_add)
        p[i] = stats.pearsonr(Y1, Y2).pvalue    # 無相関検定のp値、両側検定
        count += 1

## 結果
(p < 0.05).mean()

【実行結果】

有意水準$${\alpha=0.05}$$に設定しておきながら、事後的なデータ水増しを繰り返した結果、有意になる確率が 18% を超えました。
有意水準もマシマシですね。。。

【コード補足】np.append() と numpy ndarray
データ水増しに用いている np.append() は numpy配列 ndarray の末尾に要素を追加するコードです。
numpy配列 ndarray は numpy ライブラリで提供している配列データの形式です。
Y1 = np.append(Y1, Y1_add) で numpy配列 Y1 の末尾に Y1_add を追加して、Y1 に入れて直しています。
カッコ内の Y1 が numpy配列になっているのは、
　Y1 = rng.normal(size=n, loc=0, scale=1)
で Y1 に正規分布乱数を入れる際に、Y1 を numpy配列にしているからです。

今回の写経は以上です。
シミュレーションで統計の指標や数値と身近になれそうな予感がいたします！

シリーズの記事

次の記事

前の記事

ブログの紹介

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

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