1 確率変数.

検定や推定の理解の上で、もっとも重要なのは確率変数の理解です。
調査によって集めたデーターは確率変数ではありません。
それは、確率変数が出力した値です。
確率変数とは、データーを生み出す何かなのです。

あなた自身は確率変数です。
あなたが何か試行をするとき、ひとつのデーターが生み出されます。
同じ状況、同じ条件ならば、結果はいつも同じでしょうか？

那須の与一
あなたは弓の名手です。
的の一点に向かって矢を射込みます。
いつも的のど真ん中に命中するでしょうか？

あなたが那須の与一のような弓の名手であったとしても、矢がいつも同じポイントに当たるとは限りません。
しかし、矢はそのポイントの近くに集中します。
その集中の仕方にその人の特徴が現れます。
その散らばり具合を分布といいます。

次の図はシュミュレーションです。 扇の的ではなく数直線を的に、原点を狙って５０回矢を射込みました。

ど真ん中に命中とまではいかなくても、ほぼその近くに集まっていることがわかります。
０を中心に左右一寸（いっすん）の幅に集まり、３寸（３すん）はずれることはないといえます。
これをヒストグラムにしたものが、次の図です。

矢が、原点とその左右に集中している様子がよく分かります。

この図を、面積が確率を表すように修正してみます。

無限回の試行を想定したとき、この図がある関数式であらわされることを発見したのがガウスでした。

その関数のグラフは、

このグラフで表される分布を正規分布といいます。

この関数は、
　　　\(\phi \left( z\right) =\dfrac{1}{\sqrt{2\pi }}\exp \left( -\dfrac{z^{2}}{2}\right)\) 　　で表されます。

このグラフは、平均が０で標準偏差が１のもので、標準正規分布といいます。

2 正規乱数

Rには正規分布に従う正規乱数を出力する関数 rnorm があります。

　　　rnorm(生成する乱数の数, 平均値, 標準偏差)

次のコードは、平均 mean=2 、標準偏差 sd=1 の正規乱数を1000個生成し、ヒストグラムを出力します。
平均を中心とし左右に 2 sd の範囲に 95% のデータが集まっていることがわかります。
（正確には 1.96 倍です）

x <- rnorm(1000,mean=2,sd=1)
hist(x)

  • 面積が確率を表すようにするには、prob=TRUE を付け加えます。
    縦軸は、Density（確率密度）となります。

  • dnorm(x, 平均値, 標準偏差) は正規分布の確率密度を出力する関数です。
    dnorm(x,2,1) は、平均 mean=2 、標準偏差 sd=1 の正規分布を表す確率密度関数です。

  • curve は、関数のグラフを描画する関数です。
    add=True でグラフは重ね書きされます。

col で色を指定できます。次の図の赤い曲線は、理論上の正規分布のグラフです。

hist(x,prob=T,col="lightblue")
curve(dnorm(x,2,1),add=TRUE,col="red")

＜まとめ>.

  • rnorm（乱数を出力）.
  • hist（ヒストグラムの描画）.
  • curve（理論曲線の描画）.

これが確率分布をシミュレーションで学ぶ基本的な流れです。

確率変数 X が 平均 μ 標準偏差 σ の正規分布に従うとき、

　　　\(X\sim N(\mu,\sigma^2)\) 　　　と書きます。

2.1 和の確率分布

\(x1\sim N(2,1^2)\)、\(x2\sim N(4,1^2)\)のとき、
　\(x3=x1+x2\) の期待値と分散はどうなるでしょうか。
まず、正規乱数でx1,x2,x3 の確率密度ヒストグラムを図示します。
どのようなことが読み取れるでしょうか？

x1 = rnorm(1000,mean=2,sd=1)
x2 = rnorm(1000,mean=4,sd=1)
x3 = x1+x2
hist(x1,prob=T,col="pink",xlim=c(-2,12))
hist(x2,prob=T,col="lightblue",add=T)
hist(x3,prob=T,col="lightgreen",add=T)

関数を用いて標準偏差と平均を求めてみましょう。
mean <- c(mean(x1),mean(x2),mean(x3)) とおくと、３つの平均をまとめて表示することができます。
c( ) は、combine（結合する）関数で、数値をまとめてベクトル表示することができます。

mean <- c(mean(x1),mean(x2),mean(x3))
mean
## [1] 2.047008 3.974216 6.021224

標準偏差をまとめて表示します。

sd <- c(sd(x1),sd(x2),sd(x3))
sd
## [1] 1.0103940 0.9665289 1.3856316

Note.
期待値については、

  • \(E(x_1+x_2)=E(x_1)+E(x_2)\)

分散については、標準偏差 sd=1 は 分散 1 を意味し、

  • \(V(x_1+x_2)=V(x_1)+V(x_2)=2\)
    \(\therefore~~sd(x_1+x_2)=\sqrt{2}\)

となっていることが読み取れます。

2.2 課題

次の確率分布の期待値と分散を、正規乱数を用いて調べましょう。

  • 正規分布の差

  • 正規分布の定数倍

2.3 まとめ

和については、

  • \(V(x_1+x_2)=V(x_1)+V(x_2)\)
  • \(E(x_1+x_2)=E(x_1)+E(x_2)\)

差については、

  • \(V(x_1-x_2)=V(x_1)+V(x_2)\)
  • \(E(x_1-x_2)=E(x_1)-E(x_2)\)

また、定数\(k\) に対し、

  • \(V(kx_1)=k^2V(x_1)\) 、\(E(kx_1)=kE(x_1)\)

が成立します。

注意.

  • \(V(x_1+x_2)=V(x_1)+V(x_2)\)
  • \(V(x_1+x_1)=V(2x_1)=4V(x_1)\)

独立な確率変数は変動が相殺され、分散は元の分散の合計になります。
同じ確率変数の和は、変動が共振して大きくなり分散は４倍になります。

このように、データの変動を実感するうえでも、正規乱数は有効です。

2.4 オプション.

ggplotパッケージを使うと、ヒストグラムを透明にし、見やすくすることができます。基本の方法に慣れてから、サイトを参考に行ってください。.

data <- data.frame(x1 = rnorm(1000,mean=2,sd=1),
                   x2 = rnorm(1000,mean=4,sd=1))
data <- mutate(data, x3=x1+x2)
df <- melt(data) 
g <- ggplot(df, aes(x = value, fill = variable))
g <- g + geom_histogram(position = "identity", alpha = 0.8)
g <- g + scale_fill_npg()
plot(g)

3 カイ二乗分布

正規乱数の二乗が従う分布について調べます。

  • 標準正規乱数 \(Z\) の二乗の分布.
  • 独立な標準正規乱数 \(Z_1,Z_2,\cdots ,Z_n\) の平方和の分布.

この分布は、正規乱数の分散の推定に関連する重要なものです。

［定義］.
n 個の独立な標準正規乱数の平方和の分布を自由度 nのカイ二乗分布といいます。

  • \(Z^2 \sim \chi^2(1)\)　（自由度１のカイ二乗分布）.

  • \(Z_1^2+Z_2^2+\cdots +Z_n^2 \sim \chi^2(n)\)　（自由度nのカイ二乗分布）

［分散との関連］.
平均0,分散1である母集団から独立に得られたデータ \(Z_1,Z_2,\cdots ,Z_n\) について、
　　　\(\dfrac{\chi^2(n)}{n}=\dfrac{Z_1^2+Z_2^2+\cdots +Z_n^2}{n}\)　は n個の実現値から求めた分散であり、
その期待値は 1 になります。
また、\(\chi^2(n)\) 自身の期待値は n になります。

3.1 正規乱数の二乗.

1000個の標準正規乱数を z 、出力された値の二乗を s とします。
zの分布と s の分布を比較するため、

　　　x軸の範囲 xlim　　　y軸の範囲 ylim.

を指定しておきます。

z <- rnorm(1000,0,1)
s <- z^2
par(mfrow=c(1,2))
hist(z,prob=T,xlim=c(-10,10),ylim=c(0,1))
hist(s,prob=T,xlim=c(-10,10),ylim=c(0,1))

mean(s)
## [1] 1.053591
  • zは、1000個の標準正規乱数（平均0, 標準偏差1 の乱数）データ.
  • sは、zの値を二乗した1000個のデータ.

s については、

  • 二乗することで負のデータは右半分に移動し、グラフの高さは約２倍になり.
  • 大きい値はより大きく、小さい値はより小さくなります.
  • mean(s) は \(z^2\) の期待値であり、z の分散に近いことが分かります。

3.2 4個の正規乱数の平方和.

4個の正規乱数は独立とします。
rnorm は独立に値を出力するので、同じ関数式でも、出力される z1, z2, z3, z4 の値は異なります。

z1 <- rnorm(1000,0,1)
z2 <- rnorm(1000,0,1)
z3 <- rnorm(1000,0,1)
z4 <- rnorm(1000,0,1)
s <- z1^2+z2^2+z3^2+z4^2
par(mfrow=c(1,2))
hist(z1^2,prob=T,xlim=c(-10,10),ylim=c(0,1))
hist(s,prob=T,xlim=c(-10,10),ylim=c(0,1))

mean(s)
## [1] 4.030432

z1, z2, z3, z4 の値が独立して動くので、右に広がった分布になります。

z1, z2, z3, z4 の分散は１で、互いに独立であることから、s の期待値は 4（分散1の4倍）になることが分かります。

個数をいきなり4個にしたのは、グラフの特徴が現れやすいからです。

3.3 n個の正規乱数の平方和.

個数が増えた場合、
　　　for ( i in 1:n ) { 　　　 }
という繰り返し構文が威力を発揮します。
{　}の中には繰り返す内容が入ります。

3.3.1 n=5 の場合

s = 0
for (i in 1:5) {
  x <- rnorm(1000,0,1)
  s <- s + x^2
}
hist(s,prob=T)
curve(dchisq(x,5),add=T,col="red")

mean(s)
## [1] 5.027802

独立な標準正規分布の平方和の従う確率密度関数を、ドイツのヘルマート（Helmert 1843-1917）が発見し、ピアソンはカイ二乗分布として適合度検定に利用しています。

  • 赤の曲線は ５個の平方和が従う確率密度関数 dchisq(x,5).
  • ５個の平方和の期待値は 5（分散１の５倍）.

一般には.

  • n個の平方和が従う確率密度関数は、dchisq(x,n).
  • n個の平方和の期待値は、n（分散１のn倍）.

ということになります。

ヘルマート行列、ヘルマート変換（測地学で重要）については、大学で学んでください。測地学は地球上の座標を設定する学問で、GPSで活用されています。

参考サイト：正規母集団の分散とカイ二乗分布.

3.4 課題

n=10 の場合（10個の正規乱数の平方和）について、

  • for 構文を用いて1000個の乱数を生成し、
  • ヒストグラムと確率密度関数を図示、
  • 10個の正規乱数の平方和の期待値を求めなさい.

3.5 まとめ.

これまでは、平均 0, 標準偏差 1 の標準正規乱数の平方和について調べました。

一般の正規乱数 \(X\sim N(\mu,\sigma^2)\) については、標準化すれば上記の理論が当てはまるので、

  • 標準正規乱数 \(\dfrac{X-\mu}{\sigma}\) の二乗の期待値は 1.
  • 独立な標準正規乱数 \(\dfrac{X_1-\mu}{\sigma},\dfrac{X_2-\mu}{\sigma},\cdots ,\dfrac{X_n-\mu}{\sigma}\) の平方和の期待値は n. 

分母を払って考えれば、

  • 正規乱数 \((X-\mu)\) の二乗の期待値は \(\sigma^2\).

  • 独立な正規乱数 \((X_1-\mu),(X_2-\mu),\cdots ,(X_n-\mu)\) の平方和の期待値は \(n~\sigma^2\).

という形になります。

これは分散の定義が.
　　　\(V(X)=E((X-\mu)^2)\) であること、

独立な確率変数の和について、.
　　　\(\begin{split} V(X_1+\cdots +X_n)&=V(X_1)+\cdots +V(X_n)\\&=E((X_1-\mu)^2)+\cdots+E((X_n-\mu)^2)\end{split}\)
　　　 が成立することに対応しています。.

3.5.1 注意.

正規乱数という確率モデルのもとでは平均と分散が既知であるので標準化は容易でした。
しかし、実際のデータ分析では平均と分散は分かっておらず、標準化は不可能です。
無理に標本平均と標本分散を用いて標準化を行えば、上記の公式とはズレが生じるはずです。
「このズレを補正し分析をできるだけ正確にしよう」というのが、次章で扱う推測統計です。