正規分布

広義の定義

データが平均値の付近に集積するような分布を表します。

正規分布は、統計理論の原理を支える重要な概念です。

一般に、正規分布と表現されているものは、度数分布や確率密度関数などを対象としていることが多いです。

一定のルールやパターン（規則性）を持つデータは、十分な量の観測を行えば、度数分布が正規分布になります。

例えば、運動会の玉入れを何千〜何万回と繰り返して、スコアの度数分布を描くと正規分布になります。

これは中心極限定理といい、統計学において重要な原理原則です。

この原理を利用すれば、身の回りで起こっていることが、ものすごく小さな確率で偶然起こったものなのか、何かのパターンによって起こりうるものなのかを統計解析できます。

正規分布以外にもロングテール型、一様分布や二つ山分布などがあります。

以降、身長データを例に解説していきます。

　 ``` # 正規分布、ロングテール、一様分布、二つ山分布 import numpy as np import matplotlib.pyplot as plt norm_distY = [155,160,160,165,165,165,170,170,175] long_distY = [155,160,160,160,165,165,170,170,175] plateau_distY = [155,155,160,160,165,165,170,170,175,175] two_peak_distY = [155,155,160,160,160,165,170,170,170,175] #可視化 fig = plt.figure(figsize=(10,7)) #ヒストグラム plt.subplot(221) plt.title("norm_dist") plt.hist(norm_distY,bins=5,density=False) plt.ylim(0,3) plt.subplot(222) plt.title("long_dist") plt.hist(long_distY,bins=5,density=False,color="r") plt.ylim(0,3) plt.subplot(223) plt.title("plateau_dist") plt.hist(plateau_distY,bins=5,density=False,color="c") plt.ylim(0,3) plt.subplot(224) plt.title("two-peak_dist") plt.hist(two_peak_distY,bins=5,density=False,color="gray") plt.ylim(0,3) plt.tight_layout() plt.show() ```

正規分布の性質

偶然か起こりうるかを分ける

左右対称の綺麗な正規分布は、平均値（μ）と標準偏差（σ:sigma）から、「偶然のもの」と、「起こりうるもの」を分類するための指標を提供します。

これは、正規分布が、平均値と標準偏差さえ決まってしえば分布が定まる性質を持っているからです。

ここでは例として、平均165cm、標準偏差10の10000個の身長データで説明していきます。

分布の定義は、以下のように決められます。

ここで、 P(a <= X <= b)は変数Xが a からb の値を取る確率を表します。

• P(平均 - SD <= X <= 平均 + SD)  のとき、曲線下面積は0.683
• P(平均 - 2SD <= X <= 平均+2SD)  のとき、曲線下面積は0.954
• P(平均-3SD <= X <= 平均+3SD)  のとき、曲線下面積は0.997

これを図で説明すると、このようになります。

``` # 正規分布の性質1 from scipy.stats import gaussian_kde # カーネル密度推定（確率密度関数をガウシアン関数で推定された確率密度を作成） import numpy as np import matplotlib.pyplot as plt np.random.seed(123) mu = 165 sigma = 10 # 身長のデータセット x = np.random.normal( loc = mu, # 平均 scale = sigma, # 標準偏差 size = 10000,# 出力配列のサイズ(タプルも可) ) x = np.round(x,1) print("取得データセット数",len(x)) mean = np.mean(x) sigma = np.std(x) print("real-mean",mean,"sigma",sigma) # ヒストグラム hist, bins = np.histogram(x, bins=1000) # カーネル密度関数を作成するためのkdeオブジェクトを初期化 kde = gaussian_kde(x) # ヒストグラムの度数ではなくデータを渡す # 標本点作成 samples = np.linspace(np.min(x), np.max(x), num=len(hist)) # 密度推定 estimates = kde(samples) # 68.3%信頼区間 min_ci_683 = mu-sigma max_ci_683 = mu+sigma # 95.4%信頼区間 min_ci_954 = mu-2*sigma max_ci_954 = mu+2*sigma # 99.3%信頼区間 min_ci_993 = mu-3*sigma max_ci_993 = mu+3*sigma # min_ciはbinsのどこにあるか pos_min_683 = len(bins[bins <= min_ci_683]) pos_max_683 = len(bins[bins <= max_ci_683]) pos_min_954 = len(bins[bins <= min_ci_954]) pos_max_954 = len(bins[bins <= max_ci_954]) pos_min_993 = len(bins[bins <= min_ci_993]) pos_max_993 = len(bins[bins <= max_ci_993]) y_ci_683 = estimates[pos_min_683:pos_max_683] x_ci_683 = bins[pos_min_683:pos_max_683] y_ci_954 = estimates[pos_min_954:pos_max_954] x_ci_954 = bins[pos_min_954:pos_max_954] y_ci_993 = estimates[pos_min_993:pos_max_993] x_ci_993 = bins[pos_min_993:pos_max_993] #可視化 fig = plt.figure(figsize=(7,12)) # 68.3% ax1 = fig.add_subplot(311) ax1.set_title("68.3% confidence interval:mu-sigma") ax1.bar(x=samples, height=hist, width=0.3, alpha=0.2) ax1.set_xlim(np.min(x)-1, np.max(x)+1) ax1.set_ylim(0, np.max(hist)) #確率密度関数（ガウシアンカーネルフィッティング） ax2 = ax1.twinx() # y軸を追加 ax2.plot(samples, estimates,label='kde',color="b")#kde(確率変数リスト)でリストで指定した確率変数の確率密度を計算 ax2.set_xlim(np.min(x)-1, np.max(x)+1) ax2.set_ylim(0,np.max(estimates)) # 密度関数グラフのy軸の最大最小を設定 ax2.fill_between(x_ci_683,y_ci_683,0, alpha=0.8, color='c')#信頼区間領域を塗りつぶす ax2.fill_between(samples, estimates,0, alpha=0.2,color="m")# すべて塗りつぶす # 95.4% ax1 = fig.add_subplot(312) ax1.set_title("95.4% confidence interval:mu-2*sigma") ax1.bar(x=samples, height=hist, width=0.3, alpha=0.2) ax1.set_xlim(np.min(x)-1, np.max(x)+1) ax1.set_ylim(0, np.max(hist)) #確率密度関数（ガウシアンカーネルフィッティング） ax2 = ax1.twinx() # y軸を追加 ax2.plot(samples, estimates, label='kde',color="b")#kde(確率変数リスト)でリストで指定した確率変数の確率密度を計算 ax2.set_xlim(np.min(x)-1, np.max(x)+1) ax2.set_ylim(0,np.max(estimates)) # 密度関数グラフのy軸の最大最小を設定 ax2.fill_between(x_ci_954,y_ci_954,0, alpha=0.8, color='c')#信頼区間領域を塗りつぶす ax2.fill_between(samples, estimates,0, alpha=0.2,color="m")# すべて塗りつぶす # 99.3% ax1 = fig.add_subplot(313) ax1.set_title("99.3% confidence interval:mu-3*sigma") ax1.bar(x=samples, height=hist, width=0.3, alpha=0.2) ax1.set_xlim(np.min(x)-1, np.max(x)+1) ax1.set_ylim(0, np.max(hist)) #確率密度関数（ガウシアンカーネルフィッティング） ax2 = ax1.twinx() # y軸を追加 ax2.plot(samples, estimates,label='kde',color="b")#kde(確率変数リスト)でリストで指定した確率変数の確率密度を計算 ax2.set_xlim(np.min(x)-1, np.max(x)+1) ax2.set_ylim(0,np.max(estimates)) # 密度関数グラフのy軸の最大最小を設定 ax2.fill_between(x_ci_993,y_ci_993,0, alpha=0.8, color='c')#信頼区間領域を塗りつぶす ax2.fill_between(samples, estimates,0, alpha=0.2,color="m")# すべて塗りつぶす plt.show() ```

本来は、完全に左右対称の正規分布にこのルールを適用しますが、そのような実測データは皆無に等しいため、あえてリアルに見える正規分布のまま説明しています。

実際の統計解析では、平均と標準偏差から完全に左右対称の正規分布を仮定して算出する場合が多いです。

この身長データは、平均 165 cm、標準偏差 10 の正規分布になっていました。
一番上の式[平均 ± SD]を利用すると、155（=165-10）cmから 175（=165+10）cmまでの間に約 68.3%の人がいることになります。

その次の式[平均 ± 2SD]では、145（=165-20）cmから 185（=165+20）cmまでの間に約 95.4%の人がいることになります。

最後の次の式[平均 ± 3SD]では、135（=165-30）cmから 195（=165+30）cmまでの間に約 99.3%の人がいることになります。

論文などでよく用いられる95%信頼区間や99%信頼区間もこれらと同様に考えることができます。

• P(平均 - 1.96SD <= X <= 平均+1.96SD)  のとき、曲線下面積は0.95
• P(平均-2.58SD <= X <= 平均+2.58SD)  のとき、曲線下面積は0.99

このように考えると、確かに、データを比べるとき、このすそ野の端にいると結論付けられたデータは、比較しているデータとは同じではない（有意に差がある）と言っている意味がわかります。

この曲線下面積は永遠に1.0にはなりません。しかし、限りなく1.0に近づきます（0.9999...）。確率を計算するときは、全曲線下面積（すなわち、全確率）を便宜上、1.0として扱います。

分散との関係

正規分布は、観測データの散らばり（ばらつき）が大きくなるに連れて、度数のピークは低く、すそ野は広くなる傾向になります。

つまり、観測データの分散が大きいと、正規分布の広がりも広くなることが予想できます。

　 ``` # 正規分布の性質2 # 分散との関係 from scipy.stats import gaussian_kde # カーネル密度推定（確率密度関数をガウシアン関数で推定された確率密度を作成） import numpy as np import matplotlib.pyplot as plt np.random.seed(123) # 身長のデータセット x_s10 = np.random.normal( loc = 165, # 平均 scale = 10, # 標準偏差 size = 10000,# 出力配列のサイズ(タプルも可) ) x_s20 = np.random.normal( loc = 165, # 平均 scale = 20, # 標準偏差 size = 10000,# 出力配列のサイズ(タプルも可) ) x_s30 = np.random.normal( loc = 165, # 平均 scale = 30, # 標準偏差 size = 10000,# 出力配列のサイズ(タプルも可) ) x_s10 = np.round(x_s10,1) x_s20 = np.round(x_s20,1) x_s30 = np.round(x_s30,1) # sigma 10のヒストグラム hist_s10, _ = np.histogram(x_s10, bins=1000) #可視化 fig = plt.figure(figsize=(7,12)) # sigma 10 ax1 = fig.add_subplot(311) ax1.set_title("sigma 10") ax1.hist(x_s10, bins=len(hist_s10),density=False) ax1.set_xlim(np.min(x_s30)-1, np.max(x_s30)+1) ax1.set_ylim(0, np.max(hist_s10)) # sigma 20 ax1 = fig.add_subplot(312) ax1.set_title("sigma 20") ax1.hist(x_s20, bins=len(hist_s10),density=False) ax1.set_xlim(np.min(x_s30)-1, np.max(x_s30)+1) ax1.set_ylim(0, np.max(hist_s10)) # sigma 30 ax1 = fig.add_subplot(313) ax1.set_title("sigma 30") ax1.hist(x_s30, bins=len(hist_s10),density=False) ax1.set_xlim(np.min(x_s30)-1, np.max(x_s30)+1) ax1.set_ylim(0, np.max(hist_s10)) plt.show() ```

標準正規分布

標準正規分布は、観測データの平均を0、分散を1となるように変換してから作成した度数分布が正規分布となっている場合の分布です。

　 ``` # 標準正規分布 import matplotlib.pyplot as plt import numpy as np from scipy.stats import norm mu = 0 # 標準化済み確率変数の算術平均 variance = 1 # 標準化済み確率変数の分散 sigma = np.sqrt(variance) x = np.linspace(mu - 3.29*sigma, mu + 3.29*sigma, 1000)# 99.9%信頼区間の分布 # x = np.linspace(norm.ppf(0.001), # norm.ppf(0.999), 1000)# このようにも書けます p_dist = norm.pdf(x, mu, sigma)# 確率密度関数 plt.plot(x, p_dist) plt.show() ```

観測データの平均を0、分散を1となるように変換する方法は標準化と呼ばれ、標準化された値はZ値などと呼ばれます。

標準化を行うことにより、単位や平均値などが異なるデータ同士を単純に比較できるようになります。

上の標準正規分布の横軸を見ると、代表値の位置が0となっており、0を中心に左右に広がっていることがわかります。これが標準化の効果です。

標準化は、（観測データ - 観測データの平均値）/ 標準偏差で求まります。

例えば、身長データの平均値が165cm、標準偏差が10であるとします。
この場合、身長が170cm以上になる確率を標準正規分布から求めるには、170を標準化して（170-165/10）、この位置までの確率密度の積算を求め、全体面積(すなわち、「1.0」)から差し引けば求まります。

しかし、実際の統計解析では、疑似データのように綺麗な正規分布を描くことはほとんどありません。

最初に示したようないろいろなタイプの分布があります。

例えば、身長データを例に、2つ山分布について考えてみます。

この例では、身長150cmと170cmに偏りのあるデータになっています。

見てのとおり、度数分布も確率密度関数も正規分布になりません。

``` from scipy.stats import gaussian_kde # カーネル密度推定（確率密度関数をガウシアン関数で推定された確率密度を作成） import numpy as np import matplotlib.pyplot as plt import random # 偏りのある身長のデータセット ''' 150近辺が平均のデータ群と 170近辺が平均のデータ群を用意 ''' random.seed(123) a = [random.gauss(150, 5) for _ in range(10000)] b = [random.gauss(170, 5) for _ in range(30000)] x = a + b # リストを結合 # ヒストグラム hist, bins = np.histogram(x, bins=1000) # カーネル密度関数を作成するためのkdeオブジェクト kde = gaussian_kde(x) # 標本点を作成 pv = np.linspace(np.min(x), np.max(x), len(hist)) #可視化 fig = plt.figure() #ヒストグラム ax1 = fig.add_subplot(111) ax1.hist(x,bins=len(hist)) ax1.set_xlim(np.min(x), np.max(x)) #確率密度関数（ガウシアンカーネルフィッティング） ax2 = ax1.twinx() # y軸を追加 ax2.plot(pv,kde(pv),label='kde',color="y")#kde(確率変数リスト)でリストで指定した確率変数の確率密度を計算 ax2.set_xlim(np.min(x), np.max(x)) ax2.set_ylim(0,np.max(kde(pv))) # 密度関数グラフのy軸の最大最小を設定 plt.show() ```

標準化をしても分布の形は同じです。

``` # 続けて x_std = (x-np.mean(x))/np.std(x) # ヒストグラム hist, bins = np.histogram(x_std, bins=1000) # カーネル密度関数を作成するためのkdeオブジェクト kde = gaussian_kde(x_std) # 標本点を作成 pv = np.linspace(np.min(x_std), np.max(x_std), len(hist)) #可視化 fig = plt.figure() #ヒストグラム ax1 = fig.add_subplot(111) ax1.hist(x_std,bins=len(hist)) ax1.set_xlim(np.min(x_std), np.max(x_std)) #確率密度関数（ガウシアンカーネルフィッティング） ax2 = ax1.twinx() # y軸を追加 ax2.plot(pv,kde(pv),label='kde',color="y")#kde(確率変数リスト)でリストで指定した確率変数の確率密度を計算 ax2.set_xlim(np.min(x_std), np.max(x_std)) ax2.set_ylim(0,np.max(kde(pv))) # 密度関数グラフのy軸の最大最小を設定 plt.show() ```

正規分布とならなければ、中心極限定理を適用できません。

この度数分布のまま確率密度関数を計算し、そこから確率を算出しても、「たまたまそのデータでそうなっただけ」と一笑に付されかねません。

このようなとき、十分な観測データを収集できるまで長い年月をかけるしか無いのでしょうか。

このような課題に対する1つの対応策として、統計解析をする人たち（限られた手持ちのデータから規則性を見出そうとする人たち）は、「十分な量の観測を行えば、度数の分布が正規分布になる」という原理原則を信じて突き進むという手段を使う術を持っています。

すなわち、どのような分布であっても、一旦、正規分布を仮定するということです。

（※実際には、正規分布かどうかを検定を用いて調べてから決めます。）

どのようにするかというと、正規分布の性質である平均（代表値）と標準偏差が決まれば分布が求まる法則を利用して、どのような観測データでも正規分布を仮定した場合には標準正規分布を適用するということをやります。

標準正規分布は、確率変数の平均を0、分散を1としたときの正規分布です。

観測データを標準化すれば、必ず平均は０、分散は１（=標準偏差も1）となりますから、どのような観測データでも標準正規分布として考えられるという理屈です。

これまでに示したガウシアンフィッティングを用いた確率密度は、度数分布をそのまま反映しているために歪んだカーブを描いていました。

一方、シミュレートされる標準正規分布は常に一定の左右対称の分布になります。

標準正規分布を適用する場合、確率は標準正規分布の確率密度の曲線下面積から算出されます。

確率の算出方法は、一般の確率密度から求める場合と変わります。

一般の確率密度から求める方法（確率密度のページを参照）では、観測値が、例えば、165cmならそのまま165として、そのままのスケールで扱っていましたが、標準正規分布を使う場合は、①実測値を標準化してから、②確率（標準正規分布の確率密度から求められた面積）を計算するという手順になります。

標準化する必要がある以外は、確率の求め方は一般の確率密度から求める場合と同じです。

``` # 続けて # 標準正規分布による確率の算出 mean = np.mean(x) sigma = np.std(x) print("mean",mean,"sigma",sigma) ''' 正規分布を仮定 ''' #例えば、170センチより大きい確率は？ a = 170 za = (a-mean)/sigma # 標準化 ''' locは代表値。一般には確率変数の平均値。 scaleは標準偏差。 ※観測データを標準化した場合は、常にloc=0,scale=1 ''' pa = norm.cdf(za,loc=0,scale=1)# 標準正規分布からの累積密度（身長が170までになる確率） print("170センチより大きい確率",1-pa) # 0.30 #150～170センチである確率は？ a = 150 b = 170 # 標準化スケールへ変換 za = (a-mean)/sigma zb = (b-mean)/sigma pa = norm.cdf(za,loc=0,scale=1)# 標準正規分布からの累積密度 pb = norm.cdf(zb,loc=0,scale=1)# 標準正規分布からの累積密度 p = pb - pa print("150～170センチである確率",p) # 0.62 ```

補足

正規分布が仮定できるときに用いる統計解析手法は、パラメトリックな手法と呼ばれます。

正規分布が仮定できない場合は、平均値の信頼性が下がるために順序統計量による解析を行います。このような統計解析手法はノンパラメトリックな手法と呼ばれます。

Visionary Imaging Services, Inc.

Imaging CROサービスをご提供させていただきます。

お問い合わせお待ちしております。

https://www.vis-ionary.com/