新生児体重分布と正規分布のグラフ描画
Spiegelhalter (2019)は、”The distribution of birth weights of 1,096,277 children of non-Hispanic white women in US in 2013, born at 39-40 weeks’ gestation”のグラフを扱っている。Spiegelhalter(2019), p. xviに記載されているウェブサイトからダウンロードしたファイルを利用・参考にして、リスト１に示すPythonスクリプトを作成した。
データは、級（bin）境界値と級内の度数をスクリプトで以下のように設定している。
weights = np.array([1500, 2000, 2500,
3000, 3500, 4000, 4500, 5000])
# binの境界値
mids = weights + 250
# 中央値
n = np.array([5+48+308+1130, 12679,
124209, 442891, 389275, 108886, 14936,1345]) # bin内の度数
これらの値は、Spiegelhalter(2019), p. xviに記載されているウェブサイトからダウンロードしたRスクリプトに書かれているものである。
リスト１を実行すると、図１のグラフが表示される。
図１
左上のヒストグラムは、データが級内の度数で与えられているので、棒グラフを描く関数barによって、次のコードで描いている。
plt.bar(weights,
# 横軸
(n /
N) / 500,
# 度数を密度に変換
width=500, align='edge', edgecolor ='k', linewidth=2, alpha= 0.5)
正規分布のグラフと高さを合わせるために、度数nを次式のように密度に変換している。
(n /
N) / 500,
# 度数を密度に変換
正規分布は、データの平均値と標準偏差を正規分布の平均値と標準偏差として計算している。
すなわち、次のコードでデータに対する正規分布のグラフの描画を行っている。
x_coord = np.linspace(1500, 5500,
1000)
y_coord = ss.norm.pdf(x_coord,
birth_mean, birth_sd)
plt.plot(x_coord, y_coord, lw=3,
c='b')
なお、画面全体をk行、h列に分割して、j番目の分割画面を描画画面とすることは、次のコード
plt.subplot(khj)
で行っている。
岡本安晴（2019） いまさら聞けないPythonでデータ分析．丸善出版
Spiegelhalter, D. (2019). The art of statistics: Learning from data. Pelican
リスト１ 新生児体重の分布と正規分布の描画（圧縮ファイルbwscript.zipにまとめた）
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as ss
#
# Data
from the code downloaded from Spiegelhalter(2019)'s website
#
weights = np.array([1500, 2000, 2500,
3000, 3500, 4000, 4500, 5000])
# binの境界値
mids = weights + 250
# 中央値
n = np.array([5+48+308+1130, 12679,
124209, 442891, 389275, 108886, 14936,1345]) # bin内の度数
N = np.sum(n)
# total number
of babies
birth_mean = np.sum(mids * n) / N
# 平均値
birth_sd = (np.sum(((mids -
birth_mean)**2) * n) / N) ** 0.5
# 標準偏差
#
#
データを平均値birth_mean、標準偏差birth_sdの正規分布で表す
#
# 描画画面の準備．２×２のサブプロット画面とする
#
plt.figure(figsize=(12,7)) # 描画画面全体の大きさ
#
#
左上の画面に新生児体重のヒストグラムを表示
#
plt.subplot(221)
#
#
関数 bar でヒストグラムを描く
#
plt.bar(weights,
# 横軸
(n /
N) / 500,
# 度数を密度に変換
width=500, align='edge', edgecolor ='k', linewidth=2, alpha= 0.5)
#
#
平均birth_mean、標準偏差birth_sdの正規分布の描画
#
x_coord = np.linspace(1500, 5500,
1000)
y_coord = ss.norm.pdf(x_coord,
birth_mean, birth_sd)
plt.plot(x_coord, y_coord, lw=3,
c='b')
#
Spiegelhalter(2019)で扱われた例（p.85）
xw = 2910
plt.plot([xw,xw],
[0,ss.norm.pdf(birth_mean, birth_mean, birth_sd)],
c='r', lw=2)
#
# グラフの y=0 の位置とタイトルの表示
#
plt.yticks([0], [0])
plt.hlines(0, x_coord[0], x_coord[-1],
color='k', lw=1, alpha=0.5)
plt.title('(a) Distribution of
birthweights')
#
# 右上の画面に正規分布の平均値と±１SD、±２SD、±３SDの位置を表示
#
plt.subplot(222)
#
#
正規分布の表示
#
x_coord = np.linspace(1500, 5500,
1000)
y_coord = ss.norm.pdf(x_coord,
birth_mean, birth_sd)
plt.plot(x_coord, y_coord, lw=3,
c='b')
xw = 2910
plt.plot([xw,xw],
[0,ss.norm.pdf(birth_mean, birth_mean, birth_sd)],
c='r', lw=2)
#
# 平均値と±１SD、±２SD、±３SDの位置を垂直線で表示
#
plt.vlines([birth_mean], 0.0,
ss.norm.pdf(birth_mean, birth_mean, birth_sd),
color='k', lw=2)
plt.text(birth_mean,
ss.norm.pdf(birth_mean, birth_mean, birth_sd),
'mean', ha='center')
SDs = np.array([1, 2, 3]) * birth_sd
for i in range(3):
h = ss.norm.pdf(birth_mean +
SDs[i], birth_mean, birth_sd)
plt.vlines([birth_mean +
SDs[i], birth_mean - SDs[i]], 0.0, h,
color='k', lw=2)
plt.text(birth_mean +
SDs[i], h, f'+{i+1} SD', ha='left')
plt.text(birth_mean -
SDs[i], h, f'-{i+1} SD', ha='right')
#
# グラフの y=0 の位置とタイトルの表示
#
plt.yticks([0], [0])
plt.hlines(0, x_coord[0], x_coord[-1],
color='k', lw=1, alpha=0.5)
plt.title('(b) Mean +/- 1,2,3 SDs')
#
# 左下の画面に正規分布のパーセンタイル値の位置を表示
#
plt.subplot(223)
#
# 正規分布のグラフの表示
#
x_coord = np.linspace(1500, 5500,
1000)
y_coord = ss.norm.pdf(x_coord,
birth_mean, birth_sd)
plt.plot(x_coord, y_coord, lw=3,
c='b')
#
Spiegelhalter(2019)で扱われた例（p.85）
xw = 2910
plt.plot([xw,xw], [0,ss.norm.pdf(birth_mean,
birth_mean, birth_sd)],
c='r', lw=2)
#
#
パーセンタイル値の位置の表示
#
Pcntls = np.array([0.01, 0.05, 0.25,
0.5, 0.75, 0.95, 0.99])
x_pcntls = ss.norm.ppf(Pcntls,
birth_mean, birth_sd)
# パーセンタイル値
y_pcntls = ss.norm.pdf(x_pcntls,
birth_mean, birth_sd)
# パーセンタイル値における正規分布密度関数値
for i in range(len(x_pcntls)):
plt.vlines(x_pcntls[i], 0.0,
y_pcntls[i],
linestyles='--', color='k', lw=2)
plt.text(x_pcntls[i], 0,
f'{round(Pcntls[i]*100)}%',
va='top',
ha='center')
#
# グラフの y=0 の位置とタイトルの表示
#
plt.yticks([0], [0])
plt.hlines(0, x_coord[0], x_coord[-1],
color='k', lw=1, alpha=0.5)
plt.title('(c) Percentiles')
#
#
X軸のラベル表示
#
plt.xlabel('Birthweight (grms)',
fontsize=14)
#
# 右下の画面に正規分布における左側累積確率の表示２例
#
plt.subplot(224)
#
# 正規分布のグラフの表示
#
x_coord = np.linspace(1500, 5500,
1000)
y_coord = ss.norm.pdf(x_coord,
birth_mean, birth_sd)
plt.plot(x_coord, y_coord, lw=3,
c='b')
#
Spiegelhalter(2019)で扱われた例（p.85）
xw = 2910
plt.plot([xw,xw],
[0,ss.norm.pdf(birth_mean, birth_mean, birth_sd)],
c='r', lw=2)
#
#
2500g以下の確率
#
x_sub = np.linspace(1500, 2500, 100)
y_sub = ss.norm.pdf(x_sub, birth_mean,
birth_sd)
plt.fill_between(x_sub, y_sub,
color='aqua', alpha=0.3)
plt.text(2000, 0,
'Proportion\nbelow\n2500 gms' +
f'\n{ss.norm.cdf(2500, birth_mean, birth_sd)*100:.1f} %',
ha='center')
#
#
2910g以下の確率
#
x_sub1 = np.linspace(1500, 2910, 100)
y_sub1 = ss.norm.pdf(x_sub1,
birth_mean, birth_sd)
plt.fill_between(x_sub1, y_sub1,
facecolor=(0,0,0,0),
edgecolor='r', hatch='//')
plt.text(3400, 0,
'Proportion\nbelow\n2910 gms' +
f'\n{ss.norm.cdf(2910, birth_mean, birth_sd)*100:.0f} %',
ha='center')
#
# グラフの y=0 の位置とタイトルの表示
#
plt.yticks([0], [0])
plt.hlines(0, x_coord[0], x_coord[-1],
color='k', lw=1, alpha=0.5)
plt.title('(d) Low birth weight')
#
#
X軸のラベル表示
#
plt.xlabel('Birthweight (grms)',
fontsize=14)
plt.tight_layout()
plt.show()