「数値シミュレーションで読み解く統計のしくみ」をPythonで写経 ～ 第３章3.2「確率変数の期待値と分散」
第３章「乱数生成シミュレーションの基礎」
書籍の著者 小杉考司 先生、紀ノ定保礼 先生、清水裕士 先生
この記事は、テキスト「数値シミュレーションで読み解く統計のしくみ」第３章「乱数生成シミュレーションの基礎」3.2節「確率変数の期待値と分散」の Python写経活動 を取り扱います。
前回から引き続き、シミュレーションの核である乱数生成に取り組みます。
今回は確率分布との親交（理解）を深めるために、いくつかの確率分布に従う確率変数の期待値と分散をシミュレーションなどを使って確認します。
R の基本的なプログラム記法はざっくり Python の記法と近いです。
コードの文字の細部をなぞって、R と Python の両方に接近してみましょう！
ではテキストを開いて準備体操の旅に出発です🚀
はじめに
テキスト「数値シミュレーションで読み解く統計のしくみ」のご紹介
このシリーズは書籍「数値シミュレーションで読み解く統計のしくみ〜Rでためしてわかる心理統計」（技術評論社、「テキスト」と呼びます）の Python写経です。
テキストは、2023年9月に発売され、副題「Rでためしてわかる心理統計」のとおり、統計処理に定評のある R の具体的なコードを用いて、心理統計の理解に役立つ数値シミュレーションを実践する素晴らしい書籍です。
引用表記
この記事は、出典に記載の書籍に掲載された文章及びコードを引用し、適宜、掲載文章とコードを改変して書いています。
【出典】
「数値シミュレーションで読み解く統計のしくみ〜Rでためしてわかる心理統計」初版第１刷、著者 小杉考司・紀ノ定保礼・清水裕士、技術評論社
記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！
3.2 確率変数の期待値と分散
テキストのとおり、確率変数の期待値・分散、正規分布の再生性に取り組みます。
この記事は Jupyter Notebook 形式（拡張子 .ipynb）でPythonコードを書きます。
概ね確率分布の特性値算出には scipy.stats を用い、乱数生成には numpy.random.generator を用いるようにしています。
主に利用するライブラリをインポートします。
### インポート
# 数値計算
import numpy as np
# 確率・統計
import scipy.stats as stats
# 数値積分
from scipy import integrate
# 描画
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'Meiryo'
3.2.1 期待値
■ イカサマサイコロの期待値
3.1.1項で出現したイカサマサイコロの期待値を算出します。
１～４の出る確率がそれぞれ 0.1、５と６の出る確率がそれぞれ 0.3 です。
### 83ページ 3.1.1項のイカサマサイコロの期待値
(0.1 * 1) + (0.1 * 2) + (0.1 * 3) + (0.1 * 4) + (0.3 * 5) + (0.3 * 6)【実行結果】
4.3 です。
■ 均等なサイコロの期待値
出る目が均等に 1/6 のサイコロの期待値は、１～６の範囲の離散一様分布に従う確率変数の期待値です。
### 84ページ 公正なサイコロの期待値 a=1, b=6の離散一様分布に従う
# 設定
low = 1 # 下限
high = 6 # 上限
x = np.arange(1, 7) # 出目
# 期待値の計算
sum(stats.randint.pmf(k=x, low=low, high=high + 1) * x)【実行結果】
3.5 です。
■ 連続一様分布に従う確率変数の期待値
連続型確率分布に従う確率変数の期待値は積分計算が必要そうです。
下端$${\alpha=1}$$, 上端$${\beta=6}$$の連続一様分布に関する数値積分版の期待値算出コードです。
### 84ページ 連続一様分布に従う確率変数の期待値 数値積分版
def d_unif_exp(x, alpha=1, beta=6):
return stats.uniform.pdf(x=x, loc=alpha, scale=beta - alpha) * x
integrate.quad(func=d_unif_exp, a=1, b=6) # 数値積分の実行(積分値, 誤差推定値)【実行結果】
左が期待値です。期待値は 3.5 です。
scipy.stats の機能で連続一様分布に従う確率変数の期待値を算出します。
連続一様分布 uniform に対して mean を実行します。
### scipy.statsで期待値を算出
# パラメータ設定
alpha = 1
beta = 6
# 期待値算出
stats.uniform.mean(loc=alpha, scale=beta - alpha)【実行結果】
■ 標準正規分布に従う確率変数の期待値
テキストの実装を踏襲して期待値を計算します。
### 84ページ 標準正規分布に従う確率変数の期待値=平均パラメータ
def std_norm_exp(x, mu=0, sigma=1):
return stats.norm.pdf(x=x, loc=mu, scale=sigma) * x
integrate.quad(func=std_norm_exp, a=-np.inf, b=np.inf)【実行結果】
期待値は 0 です。
標準正規分布は平均パラメータ ０ なので、期待値は０です。
シミュレーション発動です！
標準正規分布乱数を大量に発生させて乱数の平均値をとることで、期待値をシミュレーションします。
### 85ページ 乱数生成シミュレーションで近似的に期待値を算出
rng = np.random.default_rng(seed=123)
rng.normal(size=100000, loc=0, scale=1).mean()【実行結果】
ほぼ０です。
scipy.stats の機能で標準正規分布に従う確率変数の期待値を算出します。
正規分布 norm に対して mean を実行します。
### scipy.statsで期待値を算出
# パラメータ設定
mu = 0
sigma = 1
# 期待値算出
stats.norm.mean(loc=mu, scale=sigma)【実行結果】
やっぱり０です。
3.2.2 分散
続いて、分散の算出です。
■ 離散一様分布に従う確率変数の分散
$${\alpha=1}$$、$${\beta=6}$$の離散一様分布の分散です。
テキストのコードにならって描きます。
### 85ページ 離散一様分布に従う確率変数の分散
def d_unif_var(x, alpha=1, beta=6):
# α = 1, β = 6の離散一様分布に従う確率変数の期待値
expected_val = np.arange(alpha, beta+1).mean()
var = stats.randint.pmf(k=x, low=alpha, high=beta+1) * (x - expected_val)**2
return var
d_unif_var(range(1, 7)).sum()【実行結果】
続いてテキストの通り、離散一様分布に従う確率変数の分散の公式を用いて分散を計算します。
### 86ページ 離散一様分布に従う確率変数の分散（公式利用）
alpha = 1
beta = 6
((beta - alpha + 1)**2 - 1) / 12【実行結果】
もちろん、scipy.stats でも分散を求められます。
### scipy.statsで分散を算出
# パラメータ設定
alpha = 1
beta = 6
# 分散算出
stats.randint.var(low=alpha, high=beta + 1)【実行結果】
■ 標準正規分布に従う確率変数の分散
テキストのコードにならって描きます。
### 86ページ 標準正規分布に従う確率変数の分散
def d_norm_var(x, mu=0, sigma=1):
expected_val = mu # 標準正規分布に従う確率変数の期待値
return stats.norm.pdf(x=x, loc=mu, scale=sigma) * (x - expected_val)**2
integrate.quad(func=d_norm_var, a=-np.inf, b=np.inf)【実行結果】
左が分散です。分散は１です。
標準正規分布の分散パラメータは１なので、分散は１です。
シミュレーション発動です！
標準正規分布乱数を大量に発生させて分散を算出するシミュレーションを行います。
### 86ページ 乱数生成シミュレーションで近似的に期待値を算出
# 2章で定義した自作関数
def var_p(x): # xはnumpy配列
n = len(x)
mean_x = x.mean()
var_x = sum((x - mean_x)**2) / n
return var_x
rng = np.random.default_rng(seed=0)
var_p(rng.normal(size=100000, loc=0, scale=1))【実行結果】
ほぼ１です。
もちろん、scipy.stats でも分散を求められます。
### scipy.statsで分散を算出
# パラメータ設定
mu = 0
sigma = 1
# 期待値算出
stats.norm.var(loc=mu, scale=sigma)【実行結果】
やっぱり１です。
コラム：正規分布の再生性
テキストにならって正規分布の再生性のシミュレーションを行います。
・確率変数$${X}$$：平均$${0}$$、分散$${10^2}$$の正規分布に従う
・確率変数$${Y}$$：平均$${5}$$、分散$${5^2}$$の正規分布に従う
・確率変数$${Z}$$：$${Z=X+Y}$$
・$${Z}$$は正規分布の再生性により、平均$${0+5=5}$$、分散$${10^2+5^2=125}$$の正規分布に従う（はず）
$${X}$$と$${Y}$$の正規分布乱数を加算して$${Z}$$を作成し、$${Z}$$のヒストグラムを描画します。
このヒストグラムはシミュレーションによる実現値に基づくものです。
このヒストグラムに理論値としての「平均$${5}$$、分散$${125}$$の正規分布の確率密度関数」を重ねて、一致するかどうかを確認します。
### 87ページ 正規分布の再生性のシミュレーション 88ページ図3.15
## パラメータの設定
mu_x = 0
sigma_x = 10
mu_y = 5
sigma_y = 5
n = 20000 # 生成する乱数の個数
## 確率変数x,y,zの算出
rng = np.random.default_rng(seed=0)
# 確率変数x μ = 0, σ = 10の正規分布に従う乱数
x = rng.normal(size=n, loc=mu_x, scale=sigma_x)
# 確率変数y μ = 5, σ = 5の正規分布に従う乱数
y = rng.normal(size=n, loc=mu_y, scale=sigma_y)
# 確率変数z
z = x + y
## 描画処理
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 3))
# zのヒストグラムの描画
ax.hist(z, bins=30, density=True, color='tab:blue', ec='white', alpha=0.7)
# zが従う正規分布の確率密度関数（理論値）の描画（赤い線）
line_x = np.linspace(-70, 70, 1001)
ax.plot(line_x, stats.norm.pdf(line_x, loc=mu_x+mu_y,
scale=np.sqrt(sigma_x**2 + sigma_y**2)),
color='tab:red', lw=2)
# 修飾
ax.set(xlim=(-70, 70), xlabel='z', ylabel='Density', title='Histogram of z')
ax.grid(lw=0.5);【実行結果】
変数$${Z}$$のヒストグラムと平均$${5}$$、分散$${125}$$の正規分布の確率密度関数（赤い線）はほぼ一致しました！
変数$${Z}$$の平均と分散を求めます。
### 88ページ zの平均（期待値）
z.mean()【実行結果】
平均の理論値５に近似しています。
### 88ページ zの分散
var_p(z) # 分散（2章で定義した自作関数を使用）【実行結果】
分散の理論値に近似しています。
ちなみにnumpyでも分散を計算できます。
### numpyで標本分散を算出
z.var() # ddofのデフォルト値は0【実行結果】
■ 同一の正規分布に独立に従う確率変数の標本平均が従う正規分布
標本平均の大切な性質です。
テキストの定理をお借りいたします。
平均$${\mu}$$、分散$${\sigma^2}$$に独立に従う確率変数$${X_1, X_2, \cdots, X_n}$$の標本平均$${\bar{X}}$$は、平均$${\mu}$$、分散$${\frac{\sigma^2}{n}}$$の正規分布に従います。
この定理をシミュレーションで確かめます。
平均$${\mu=50}$$、分散$${\sigma^2=10^2}$$の正規分布に独立に従う確率変数$${X_1, X_2, X_3, X_4}$$の標本平均$${\bar{X}}$$を10000個、乱数生成で作成し、$${\bar{X}}$$のヒストグラムを描画します。
理論値としての平均$${50}$$、分散$${10^2/4=25}$$（標準偏差$${5}$$）の正規分布の確率密度関数と重ねて、ヒストグラムとの一致を確かめます。
### 89ページ 確率変数の平均のヒストグラムの描画 90ページ図3.16
n = 4
mu = 50
sigma = 10
iter = 10000
means = np.zeros(iter)
rng = np.random.default_rng(seed=123)
for i in range(iter):
means[i] = rng.normal(size=n, loc=mu, scale=sigma).mean()
## 平均のヒストグラムの描画処理
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 3))
# 平均のヒストグラムの描画
ax.hist(means, bins=30, density=True, color='tab:blue', ec='white', alpha=0.7)
# 平均が従う正規分布の確率密度関数（理論値）の描画（赤い線）
line_x = np.linspace(means.min(), means.max(), 1001)
ax.plot(line_x, stats.norm.pdf(line_x, loc=mu, scale=sigma / np.sqrt(n)),
color='tab:red', lw=2)
# 修飾
ax.set(xlabel='mean', ylabel='Density', title='Histogram of mean')
ax.grid(lw=0.5);【実行結果】
ほぼ一致しました！
今回の写経は以上です。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
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の教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
この記事が気に入ったらサポートをしてみませんか？