2012/09/13

文系のための「正規分布」(2)

正規分布の存在意義というのは、結局のところ、平均標準偏差さえ分かれば、
釣り鐘状」に値が分布するような関数を作ることができ、その関数を用いることで、
ある値を取る確率を簡単に計算できるということであった。

もちろん、釣り鐘状には分布しないこともあって、
その時には、残念ながら正規分布を前提にすることはできないのであるが、
別の様々な分布を前提とすることになる。

他の様々な分布を知ることは重要ではあるが、ここでは正規分布に限った話にする。
まずは、正規分布の意味が解っていれば、他の分布関数を理解する上でも、
多少の手助けにはなるだろう。

さて。前回までの話では、一般的な正規分布の話をしてきたが、
今回の話では、少し特殊な場合での正規分布について整理する。

どのように特殊なのかというと、つまり、
平均が「0」で標準偏差が「1」となるような正規分布である。
この種の正規分布のことを、一般的には、「標準正規分布」と呼ぶ。

この標準正規分布を用いると、色々と計算が便利になるので、様々なとこに出てくる。
なぜ、標準正規分布を用いるか、という話はとりあえず置いておいて、まずは基本的な話から。

さて、前回の話では、一般的な正規分布の関数の式を以下のように表した。



ここで、μは平均を表し、σは標準偏差を表したのであった。
標準正規分布では、μ=0、σ=1 となるので、以下のように表される。



確か、正規分布は、平均標準偏差を与えれば良いだけなので、
これだけで、釣り鐘状の山を描くことができる。Rでは以下のようにする。
# Xの取り得る範囲を決める。とりあえず、-4~4くらいで。これはテキトー。
X.range <- seq(-4, 4, 0.01)

# 標準正規分布のX.range の間の f(x)の値を得る
X.dnorm <- dnorm(X.range)

# 山を描く
plot(X.range, X.dnorm, type="l", xlab="", ylab="", ylim=c(0,0.5))
title(main="Bell curve for Standard Normal Distribution")
title(xlab="Range of x value(Standard deviation)", ylab="Density")
abline(h=0)

# 平均を通る線を引く
abline(v=0, lty=2, col="red")
text(0, 0.45, expression(mu==0), cex=1.8)

この図において、x軸の値が標準偏差であり、y軸の値がその時の確率の密度を表す。
標準偏差「0」の所で山の頂点があり、ここが平均を示している。
また、±4」付近に山裾が存在しているのが確認できる。

標準正規分布において重要なことは、
平均からの距離を標準偏差の倍数で確率を推測するということであり、
区切りの良い距離として±1σ倍±2σ倍±3σ倍を用いることが多い。

とりあえず、この状況を可視化して整理してみる。まずは、±1σ倍の場合から。

なお、今回から、x軸の軸ラベルには、一般的な表記を用いるが、
平均「0」、標準偏差「1」の場合には、先ほどの図と同じ状況になる。
# 山を描く
plot(X.range, X.dnorm, type="l", xlab="", ylab="", xaxt="n", ylim=c(0,0.5))
axis(1, at=seq(-3, 3, by=1), labels=c("μ-3σ","μ-2σ", "μ-1σ", "μ", "μ+1σ", "μ+2σ", "μ+3σ"))
title(main="Bell curve for Standard Normal Distribution")
title(xlab="Range of x value(Standard deviation)", ylab="Density")
abline(h=0)

# 面積の部分を塗りつぶす
poly.x <- seq(-1,1,0.01)
poly.y <- c(0, dnorm(poly.x), 0)
polygon(c(-1, poly.x, 1), poly.y, col="lightblue")

# 平均を通る線を引く
abline(v=0, lty=2, col="red")

# キャプションを入れる。
prob <- paste(round(pnorm(1) - pnorm(-1), 3)*100, "%")
text(0, 0.1, prob, cex=1.8)

次は、±2σ倍の場合。

# 山を描く
plot(X.range, X.dnorm, type="l", xlab="", ylab="", xaxt="n", ylim=c(0,0.5))
axis(1, at=seq(-3, 3, by=1), labels=c("μ-3σ","μ-2σ", "μ-1σ", "μ", "μ+1σ", "μ+2σ", "μ+3σ"))
title(main="Bell curve for Standard Normal Distribution")
title(xlab="Range of x value(Standard deviation)", ylab="Density")
abline(h=0)

# 面積の部分を塗りつぶす
poly.x <- seq(-2,2,0.01)
poly.y <- c(0, dnorm(poly.x), 0)
polygon(c(-2, poly.x, 2), poly.y, col="lightblue")

# 平均を通る線を引く
abline(v=0, lty=2, col="red")

# キャプションを入れる。
prob <- paste(round(pnorm(2) - pnorm(-2), 3)*100, "%")
text(0, 0.1, prob, cex=1.8)

次は、±3σ倍の場合。

# 山を描く
plot(X.range, X.dnorm, type="l", xlab="", ylab="", xaxt="n", ylim=c(0,0.5))
axis(1, at=seq(-3, 3, by=1), labels=c("μ-3σ","μ-2σ", "μ-1σ", "μ", "μ+1σ", "μ+2σ", "μ+3σ"))
title(main="Bell curve for Standard Normal Distribution")
title(xlab="Range of x value(Standard deviation)", ylab="Density")
abline(h=0)

# 面積の部分を塗りつぶす
poly.x <- seq(-3, 3, 0.01)
poly.y <- c(0, dnorm(poly.x), 0)
polygon(c(-3, poly.x, 3), poly.y, col="lightblue")

# 平均を通る線を引く
abline(v=0, lty=2, col="red")

# キャプションを入れる。
prob <- paste(round(pnorm(3) - pnorm(-3), 3)*100, "%")
text(0, 0.1, prob, cex=1.8)

さて、ところで、平均「0」、標準偏差「1」、という状況から連想することはないだろうか?
確かにどこかで聞いたような。そうそう、標準化の話に似ている。
標準化の話では、標準偏差ではなく、分散が「1」と説明したが、結局は同じこと。

念のために。標準偏差を「」としたとき、分散は「」である。つまり、二乗ということ。
したがって、「」の場合、「」となる。

つまり、何が言いたいかというと、正規分布に従うデータがあって、
そのデータを標準化したデータは、標準正規分布に従うのである。
ここで、標準化の話を思い出してみる。確か、以下のような式であった。



重要なのは、右端の標準化の式である。
これを使えば、一般的な正規分布標準正規分布に変換できる。

なぜ、これが便利なのか?

前回の話の最期に、数学嫌いの人にとって目を背けたくなる積分の式が登場した。
これを実際に「手計算」で解くのは、少々厄介なのである。
もっとも、現在は、コンピュータが計算してくれるので問題は無いのであるが...。

そういった状況があって、コンピュータが存在していない時には、
標準正規分布表」というものを使っていた。

あらゆる正規分布の表を作ろうとすると、無限個の表を作らないといけないので、
要するに、標準正規分布の場合だけの表を作り、それを参照していたのである。
っということで、試しに標準正規分布表を作ってみる。

かなり、長くなってしまうので、マイナス方向は省略平均を中心に対称なので構わない。
単なる表を作るだけの作業なので、余裕の無い人はパスして構わない。

# 行数を数えるための変数を初期化する
nrow <- 0

# 結果を格納するための行列を準備する
N <- matrix(ncol=10, nrow=41)

# 標準正規分布表を作成する
for(i in seq(0, 4, by=0.1)){
    nrow <- nrow + 1    # 行数を一行進める
    ncol <- 0               # 列数を数えるための変数を初期化する
    for(j in seq(0, 0.09, by=0.01)){
        ncol <- ncol + 1  # 列数を一行進める

        # 平均から上の部分の面積を計算する。これが確立。
        N[nrow, ncol] <- round(pnorm(i+j) -pnorm(0), 4)
    }
}

# 列名と行名を入れる

rownames(N)<-seq(0, 4, by=0.1)
colnames(N)<-seq(0, 0.09, by=0.01)

そして、これを実行すると、標準正規分布表が出来上がる。

> N
         0   0.01   0.02   0.03   0.04   0.05   0.06   0.07   0.08   0.09
0   0.0000 0.0040 0.0080 0.0120 0.0160 0.0199 0.0239 0.0279 0.0319 0.0359
0.1 0.0398 0.0438 0.0478 0.0517 0.0557 0.0596 0.0636 0.0675 0.0714 0.0753
0.2 0.0793 0.0832 0.0871 0.0910 0.0948 0.0987 0.1026 0.1064 0.1103 0.1141
0.3 0.1179 0.1217 0.1255 0.1293 0.1331 0.1368 0.1406 0.1443 0.1480 0.1517
0.4 0.1554 0.1591 0.1628 0.1664 0.1700 0.1736 0.1772 0.1808 0.1844 0.1879
0.5 0.1915 0.1950 0.1985 0.2019 0.2054 0.2088 0.2123 0.2157 0.2190 0.2224
0.6 0.2257 0.2291 0.2324 0.2357 0.2389 0.2422 0.2454 0.2486 0.2517 0.2549
0.7 0.2580 0.2611 0.2642 0.2673 0.2704 0.2734 0.2764 0.2794 0.2823 0.2852
0.8 0.2881 0.2910 0.2939 0.2967 0.2995 0.3023 0.3051 0.3078 0.3106 0.3133
0.9 0.3159 0.3186 0.3212 0.3238 0.3264 0.3289 0.3315 0.3340 0.3365 0.3389
1   0.3413 0.3438 0.3461 0.3485 0.3508 0.3531 0.3554 0.3577 0.3599 0.3621
1.1 0.3643 0.3665 0.3686 0.3708 0.3729 0.3749 0.3770 0.3790 0.3810 0.3830
1.2 0.3849 0.3869 0.3888 0.3907 0.3925 0.3944 0.3962 0.3980 0.3997 0.4015
1.3 0.4032 0.4049 0.4066 0.4082 0.4099 0.4115 0.4131 0.4147 0.4162 0.4177
1.4 0.4192 0.4207 0.4222 0.4236 0.4251 0.4265 0.4279 0.4292 0.4306 0.4319
1.5 0.4332 0.4345 0.4357 0.4370 0.4382 0.4394 0.4406 0.4418 0.4429 0.4441
1.6 0.4452 0.4463 0.4474 0.4484 0.4495 0.4505 0.4515 0.4525 0.4535 0.4545
1.7 0.4554 0.4564 0.4573 0.4582 0.4591 0.4599 0.4608 0.4616 0.4625 0.4633
1.8 0.4641 0.4649 0.4656 0.4664 0.4671 0.4678 0.4686 0.4693 0.4699 0.4706
1.9 0.4713 0.4719 0.4726 0.4732 0.4738 0.4744 0.4750 0.4756 0.4761 0.4767
2   0.4772 0.4778 0.4783 0.4788 0.4793 0.4798 0.4803 0.4808 0.4812 0.4817
2.1 0.4821 0.4826 0.4830 0.4834 0.4838 0.4842 0.4846 0.4850 0.4854 0.4857
2.2 0.4861 0.4864 0.4868 0.4871 0.4875 0.4878 0.4881 0.4884 0.4887 0.4890
2.3 0.4893 0.4896 0.4898 0.4901 0.4904 0.4906 0.4909 0.4911 0.4913 0.4916
2.4 0.4918 0.4920 0.4922 0.4925 0.4927 0.4929 0.4931 0.4932 0.4934 0.4936
2.5 0.4938 0.4940 0.4941 0.4943 0.4945 0.4946 0.4948 0.4949 0.4951 0.4952
2.6 0.4953 0.4955 0.4956 0.4957 0.4959 0.4960 0.4961 0.4962 0.4963 0.4964
2.7 0.4965 0.4966 0.4967 0.4968 0.4969 0.4970 0.4971 0.4972 0.4973 0.4974
2.8 0.4974 0.4975 0.4976 0.4977 0.4977 0.4978 0.4979 0.4979 0.4980 0.4981
2.9 0.4981 0.4982 0.4982 0.4983 0.4984 0.4984 0.4985 0.4985 0.4986 0.4986
3   0.4987 0.4987 0.4987 0.4988 0.4988 0.4989 0.4989 0.4989 0.4990 0.4990
3.1 0.4990 0.4991 0.4991 0.4991 0.4992 0.4992 0.4992 0.4992 0.4993 0.4993
3.2 0.4993 0.4993 0.4994 0.4994 0.4994 0.4994 0.4994 0.4995 0.4995 0.4995
3.3 0.4995 0.4995 0.4995 0.4996 0.4996 0.4996 0.4996 0.4996 0.4996 0.4997
3.4 0.4997 0.4997 0.4997 0.4997 0.4997 0.4997 0.4997 0.4997 0.4997 0.4998
3.5 0.4998 0.4998 0.4998 0.4998 0.4998 0.4998 0.4998 0.4998 0.4998 0.4998
3.6 0.4998 0.4998 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999
3.7 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999
3.8 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999 0.4999
3.9 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000
4   0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000

少々見難い出力となってしまったが、一行目は行名で、一列目は列名である。
ちょっと、特殊な見方になっていて、行方向は小数点第一位を示していて、
列方向は小数点第二位を示している。

つまり、0.36σ のときの面積であれば、
0.3」+「0.06」なので、4行目の7列目の値を見る。
すると、「0.1406」となっている。これが、平均「0」からの面積となる。

さて、ここで、分を見てみると、11行目の1列目の値なので「0.3413」となっている。
この値を2倍すると「0.6826」となり、±1σの値と一致する。
つまり、標準正規分布表は以下のような状況を示している。
# 片袖の標準正規分布の図を描く
x <- seq(0, 4, by=0.01)
y <- dnorm(x)
plot(x, y, type="l", xlab="", ylab="", xaxt="n")
title(main="Bell curve for Standard Normal Distribution")
title(xlab="Range of x value(Standard deviation)", ylab="Density")
abline(h=0)

# 面積の部分を描く
poly.x <- seq(0, 1.8, 0.01)
poly.y <- c(0, dnorm(poly.x), 0)
polygon(c(0, poly.x, 1.8), poly.y, col="lightgray")
text(0.8, 0.1, expression(integral(f(x)*dx,0,a)), cex=1.8)

現在では、すっかり、標準正規分布表を用いることは無くなったが、
かつては、この表を用いて確率を求めていたのである。

さて、重要なことは、要するに正規分布に従うあらゆるデータは、
そのデータを標準化することで標準正規分布に変換して考えることができる。

0 件のコメント:

コメントを投稿