男女出生率のベイズ分析
ラプラスの男女の出生数データが、B. Lambert (2018), “A student’s guide to Bayesian Statistics”, SAGE Publications Ltdのp. 16に、女児241,945人、男児251,527人として載っており、ラプラスはこのデータから「女子の方が男子より産まれやすい」という確率が約
であると計算したと説明されている。このことを確認するPythonスクリプトを用意してみた。
確率分布として２項分布を設定する。２項分布を用いたベイズ分析については、＜岡本安晴「いまさら聞けないPythonでデータ分析」丸善出版＞、pp. 169-176などを参照されたい。
女児が産まれる確率をpとおくと、データ「女児241,945人、男児251,527人」に対する事後分布は、事前分布として一様分布を設定したとき
で与えられる。すなわち、ベータ分布
である。この分布を描くPythonスクリプトをリスト１のように用意した。
リスト１
import scipy.stats as ss
import numpy as np
import matplotlib.pyplot as plt
dist = ss.beta(241945+1, 251527+1)
x = np.linspace(0, 1, 1000)
y = dist.pdf(x)
plt.rcParams['font.family'] = ['Yu
Gothic', 'Noto Sans CJK JP', 'Hiragino Maru Gothic Pro']
plt.title('女子の出生確率pの分布', fontsize = 20)
plt.plot(x, y)
plt.vlines(0.5, 0, 500, color = 'r',
label = 'p = 0.5')
plt.xlabel('p', fontsize = 16)
plt.yticks([])
plt.legend()
plt.savefig('Fig_1.png')
plt.show()
リスト１のスクリプトを実行すると、図１のグラフが表示される。
図１
事後分布がp=0.5の近くに集まっているので、横軸を拡大して表示するため、リスト２のスクリプトを用意した。
リスト２
import scipy.stats as ss
import numpy as np
import matplotlib.pyplot as plt
dist = ss.beta(241945+1, 251527+1)
x_1 = np.linspace(0.48, 0.52, 1000)
y_1 = dist.pdf(x_1)
plt.rcParams['font.family'] = ['Yu
Gothic', 'Noto Sans CJK JP', 'Hiragino Maru Gothic Pro']
plt.title('女子の出生確率pの分布', fontsize = 20)
plt.plot(x_1, y_1)
plt.vlines(0.5, 0, 500, color = 'r',
label = 'p = 0.5')
plt.xlabel('p', fontsize = 16)
plt.yticks([])
plt.legend()
plt.savefig('Fig_2.png')
plt.show()
リスト２のスクリプトを実行すると図２のグラフを得る。
図２
事後分布がp=0.5から明確に離れていることが分かる。図２のグラフから、「p>0.5」すなわち「女児の方が産まれやすい確率」はほぼ０であると言えるが、事後分布における「p>0.5」の確率を求めてみる。非常に小さい（小数点以下の桁数が実数の精度より小さい）値であるので、男児の産まれる確率をpとしてベータ分布
の累積確率分布関数のp=0.5における値として求める。Pythonスクリプトをリスト３に示す。
リスト３
import scipy.stats as ss
dist_1 = ss.beta(251527+1, 241945+1)
print('p = ', dist_1.cdf(0.5))
リスト３の実行例を図３に示す。
図３
Laplaceの計算結果、約が確認できた。
ちなみに、女児の産まれる確率の95％確信区間を求めると
95% CI = [0.489, 0.492]
を得た。スクリプトをリスト４に示す。
リスト４
import numpy as np
import scipy.stats as ss
import scipy.optimize as so
class BirthRates:
def __init__(self, ng, nb):
self.MyBeta = ss.beta(ng + 1, nb + 1)
self.z_v = 0.0
def set_z_v(self, z_v):
self.z_v = z_v
# 累積確率分布 - z_v
def f(self, v):
return self.MyBeta.cdf(v) - self.z_v
# fの導関数
def df(self, v):
return self.MyBeta.pdf(v)
BR = BirthRates(241945, 251527)
f = BR.f
df = BR.df
x = np.linspace(0, 1, 1000)
y = df(x)
v_init = x[np.argmax(y)] # Newton法の初期値
BR.set_z_v(0.025)
# 累積確率が0.025になる値を求める。
v025 = so.newton(f, v_init, df)
BR.set_z_v(0.975)
# 累積確率が0.975になる値を求める
v975 = so.newton(f, v_init, df)
print('95% CI = [{0:.3f},
{1:.3f}]'.format(v025, v975))