【医療統計入門】Pythonで学ぶ正規分布・標準偏差・標準誤差①

はじめに

今回は「正規分布標準偏差・標準誤差」などの基本的な統計知識を記事にしていきたいと思います。
Jupyter-Notebookを起動させて必要なものをインポートしておきます。

import numpy as np
import pandas as pd
import scipy as sp
from scipy import stats
from numpy.random import normal,choice
import matplotlib.pyplot as plt
import matplotlib.animation as ani
import seaborn as sns
sns.set()
%precision 3
%matplotlib inline
標本抽出シミュレーション

1回の調査や実験のことを試行と呼びます。
1回の標本抽出は1回の試行とも言えます。
「標本分布」とは標本の統計量が従う確率分布のことです。
母集団からの標本抽出シミュレーションを1万回行ったとすると、1万個の標本が得られます。
標本からそれぞれ標本の平均が計算でき、標本平均が1万個できることになります。
この「1万個の標本平均の従う確率分布」が標本分布になります。

母集団の設定

今回の設定は「男女それぞれのヘモグロビン濃度」です。
男性は「平均13で標準偏差1」の正規分布/女性は「平均12で標準偏差0.8」の正規分布に従うものとしました。
母集団の数はそれぞれ10万人です。

np.random.seed(1)
man_data = normal(13,1,100000)
sns.distplot(man_data)

f:id:Medicmed:20180617111257p:plain

np.random.seed(1)
women_data = normal(12,0.8,100000)
sns.distplot(women_data,color='indianred')

f:id:Medicmed:20180617111401p:plain

男性・女性それぞれの母集団のグラフを重ねてグラフを作成します。

sns.distplot(man_data, label='man_data')
sns.distplot(women_data,color='indianred', label='women_data')
plt.legend()

f:id:Medicmed:20180617221637p:plain

母集団それぞれのデータをデータフレームの形に変形させます。

df_man = pd.DataFrame(man_data)
df_women = pd.DataFrame(women_data)
mf_df = pd.concat([df_man,df_women],axis=1)
mf_df.columns = ['man','women']
mf_df.head()

f:id:Medicmed:20180617111621p:plain

要約統計量を求める
mf_df.describe()

f:id:Medicmed:20180617111821p:plain

標本抽出シミュレーションを行う

サンプルサイズは10人とします。
標本抽出を1万回繰り返すシミュレーションを行います。
標本が1万個できるので、それぞれの標本の平均を求めます。

それぞれの標本平均を格納する変数を用意します。

man_sample_mean = np.zeros(10000)
women_sample_mean = np.zeros(10000)

seed()を用いて乱数を固定し、男性・女性をそれぞれ復元抽出で10人ランダムサンプリングするのを1万回行います。

np.random.seed(1)
for i in range(0,10000) :
man_sample = choice(man_data,10,replace=True)
women_sample = choice(women_data,10,replace=True)
man_sample_mean[i] = sp.mean(man_sample)
women_sample_mean[i] = sp.mean(women_sample)

1万個の標本データを男女別に格納したので、それぞれの「標本平均の平均値」を求めます。

標本平均の平均値を求める
Msample_mean_df = pd.DataFrame(man_sample_mean)
Fsample_mean_df = pd.DataFrame(women_sample_mean)
mf_sample_df = pd.concat([Msample_mean_df,Fsample_mean_df],axis=1)
mf_sample_df.columns = ['man_sample_mean','women_sample_mean']
mf_sample_df.describe()

f:id:Medicmed:20180617112630p:plain

男女それぞれの「標本平均の平均値」は母集団の平均(母平均)にほぼ一致していることが分かります。
しかし、男女それぞれの母集団の標準偏差に比べて「標本平均の標準偏差」はかなり小さい値になっていることが分かります。
「標本平均の平均値」は母集団の平均を推定するのに良い推定であるのに対して「標本平均の標準偏差」は母集団の標準偏差を推定するには過小評価してしまうことが分かりました。

「標本平均」のヒストグラムを描いてみる
sns.distplot(man_sample_mean, color='blue',label='man_sample_mean')
sns.distplot(women_sample_mean,color='red', label='women_sample_mean')
plt.legend()

f:id:Medicmed:20180617221812p:plain

最初にグラフにした「母集団」のグラフに比べて、それぞれ「バラツキ」=「標準偏差」が小さくなっていることがグラフからも分かります。
…と言っても良くわからないので母集団のグラフと標本平均のグラフを重ね合わせてみます。

母集団のグラフと標本平均のグラフを同時に描く
plt.figure(figsize=(10,8))
sns.distplot(man_data, label='man_data')
sns.distplot(women_data,color='indianred', label='women_data')
sns.distplot(man_sample_mean, color='blue',label='man_sample_mean')
sns.distplot(women_sample_mean,color='red', label='women_sample_mean')
plt.vlines(ymin=0,ymax=1.6,x=[sp.mean(man_data),sp.mean(women_data)],linestyles='dashed')
plt.legend(fontsize=15)

f:id:Medicmed:20180617222035p:plain

点線で示してあるのは、母集団それぞれの「平均値」です。
「標本平均の平均値」は、ほとんど母集団の平均値と一致していることがグラフから読み取れます。
一方で、「標本平均」のグラフはそれぞれ鋭いグラフになっていて「標本平均の標準偏差」は小さいことが分かります。
この「標本平均の標準偏差」のことを「標準誤差」と呼び、通常はSEM(Standard Error of Mean)と表示されます。
標準誤差の値は基本的に標準偏差よりも小さい値になることが知られています。

サンプルサイズと標本平均と母平均の関係

サンプルサイズが大きいなら標本平均は母平均に近くなることを確認します。
「10〜10010」まで100区切りで変化させたサンプルサイズを用意します。
標本平均を格納する変数を用意します。

size_array = np.arange(start=10, stop=100100, step=100)
sample_mean_array_size = np.zeros(len(size_array))
np.random.seed(1)
for i in range(0, len(size_array)) :
sample = choice(man_data, size=size_array[i], replace=True)
sample_mean_array_size[i] = sp.mean(sample)

前準備ができたので、グラフを作成します。

plt.figure(figsize=(10,8))
plt.plot(size_array, sample_mean_array_size)
plt.xlabel('sample size')
plt.ylabel('sample mean')

f:id:Medicmed:20180617114044p:plain

アニメーションを作成する

Matplotlibのanimationを使用して確認してみます。
今回はArtistAnimation()ではなくFuncAnimation()を使用します。

def animation_update(i) :
plt.plot(size_array[i], sample_mean_array[i], 'bo' )
fig = plt.figure(figsize=(8,8))
ani = ani.FuncAnimation(fig,animation_update,interval=1)
ani.save('hoge.gif', writer='imagemagick', fps=5)
plt.show()

f:id:Medicmed:20180617114410g:plain

サンプルサイズが大きくなるにつれて、母集団の平均値にほぼ一致していくのがアニメーションで分かります。

何故かサンプルサイズが1万の所でアニメーションが終わりますが、よく分かりません。
雰囲気が分かれば大丈夫です()

ぜひ参考にしてみてください。
終わり。

今回参考にした書籍はコチラ!

Pythonで学ぶあたらしい統計学の教科書

Pythonで学ぶあたらしい統計学の教科書

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

このサイトはスパムを低減するために Akismet を使っています。コメントデータの処理方法の詳細はこちらをご覧ください