第3章のPythonコード

第3章 確率論の基礎

サンプルデータ

wage.csv：男性労働者の賃金データ．

モジュールのインポート

import numpy as np
import math
import matplotlib.pyplot as plt
import pandas as pd
import scipy.stats as stats
import seaborn as sns

乱数の設定

np.random.seed(seed = 2022)

3.1.2 試行と事象

coin = ['Head', 'Tail']
np.random.choice(coin, size = 100, replace = True)

3.1.3 コイン投げのシミュレーション

math.factorial(100) / (math.factorial(50) ** 2) / (2 ** 100)
coin = [1, 0]
z = np.random.choice(coin, size = 100, replace = True)

z.sum()
S = 100000
rec = np.zeros(S)
coin = [1, 0]

for i in range(S):
  z = np.random.choice(coin, 100, replace = True)
  rec[i] = z.sum()
plt.hist(rec, bins = 20)
plt.show()
pd.DataFrame(rec).describe()

3.1.4 論理演算によるカウントの方法

2 > 1
2 > 1000
200 == 100 * 2
True + True
True * False
count = (rec == 50)
print(count)
count.sum()
count.mean()

3.2.3 独立ではない例

# 独立性の確認

S = 10000 # シミュレーション回数
X = np.random.normal(loc = 50, scale = 10, size = S) # Xを抽出
Y = np.random.normal(loc = 50, scale = 10, size = S) # Yを抽出
Z = X + Y # Zを構成
# Pr(X > 70) * Pr(Z > 100)
(X > 70).mean() * (Z > 100).mean() 
# Pr(X > 70 かつ Z > 100)
((X > 70) * (Z > 100)).mean()

3.2.4 独立性と相関係数

X = np.random.normal(loc = 50, scale = 10, size = 100000)
Y = np.random.normal(loc = 50, scale = 10, size = 100000)

np.corrcoef(X, Y)
Z = -((X - 50) ** 2) / 10
np.corrcoef(X, Z)
plt.scatter(X, Z)
plt.show()

3.3.1 分布関数

x = np.linspace(0, 100, 1000)
plt.plot(x, stats.norm.cdf(x, 50, 10))
plt.show()
print(stats.norm.cdf(60, 50, 10) - stats.norm.cdf(40, 50, 10))

3.3.2 確率密度関数

print(stats.norm.cdf(50 + 0.1, 50, 10) - stats.norm.cdf(50, 50, 10))
print(stats.norm.pdf(50, 50, 10))
print(stats.norm.pdf(80, 50, 10))
x = np.linspace(start = 0, stop = 100, num = 1000)
plt.plot(x, stats.norm.pdf(x, 50, 10))
plt.show()

3.3.7 データによる条件付期待値の推定

wagedata = pd.read_csv("wage.csv")
wagedata.head()
wagedata.describe()
sns.pairplot(wagedata)
educ12 = wagedata.loc[wagedata['educ'] == 12]
educ16 = wagedata.loc[wagedata['educ'] == 16]
educ12['wage'].mean()
educ16['wage'].mean()
educ11_less = wagedata.loc[wagedata['educ'] < 12]
educ12_more = wagedata.loc[wagedata['educ'] >= 12]
educ11_less['wage'].mean()
educ12_more['wage'].mean()

3.4.2 中心極限定理のシミュレーション

S = 10000
n = 10000
Zn = np.zeros(S)

for i in range(S):
  X = np.random.normal(loc = 50, scale = 10, size = n)
  Xbar = X.mean()
  Sn = X.var()
  Zn[i] = (n ** 0.5) * (Xbar - 50) / (Sn ** 0.5)

plt.hist(Zn, bins = 20)
plt.show()
pd.DataFrame(Zn).describe()

3.4.3 信頼区間のシミュレーション

# 信頼区間シミュレーション

S = 10000 # シミュレーション回数
n = 10000 # 標本の大きさ
rec = np.zeros(S) # 結果記録用のベクトル

for i in range(S):  # 繰り返し開始
  X = np.random.normal(loc = 50, scale = 10, size  = n)  # N(50,10^2) から標本抽出
  Xbar = X.mean()
  Sn = X.var()
  rec[i] = (Xbar - 1.96 * ((Sn / n) ** 0.5) < 50) * (50 < Xbar + 1.96 * ((Sn / n) ** 0.5))

rec.mean() # 確率を計算

3.4.4 信頼区間の導出

print(stats.norm.cdf(69.6, 50, 10) - stats.norm.cdf(30.4, 50, 10))