1 確率変数

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

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

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

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

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

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

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

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

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

その関数のグラフは、

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

この関数は、　　　\(\phi \left( z\right) =\dfrac{1}{\sqrt{2\pi }}\exp \left( -\dfrac{z^{2}}{2}\right)\)　　で表されます。
このグラフは、平均が０で標準偏差が１のもので、標準正規分布といいます。

1.1 正規乱数

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

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

面積が確率を表すようにするには、prob=TRUE を付け加えます。
縦軸は、Density（確率密度）となります。
赤は理論上の正規分布のグラフです。

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

ggplotを用いると近似曲線を表示することができます。

df <- data.frame(value = rnorm(1000, 2, 1))
head(df)

value
1.0387410
1.7309114
1.9833183
1.2326595
2.7530924
0.3288571

g <- ggplot(df, aes(x = value, y = ..density..))
g <- g + geom_histogram( alpha = 0.3)
g <- g + geom_density(aes(alpha = 0.5), show.legend = F)
g <- g + scale_fill_npg() + scale_color_npg()
plot(g)

## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

2 確率分布

2.1 和の確率分布

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

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)

## No id variables; using all as measure variables

g <- ggplot(df, aes(x = value, fill = variable))
g <- g + geom_histogram(position = "identity", alpha = 0.8)
g <- g + scale_fill_npg()
plot(g)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

関数を用いて標準偏差と平均を求めてみましょう。
標準偏差はsd 、平均はmean です。

x1 <- rnorm(1000,mean=2,sd=1)
x2 <- rnorm(1000,mean=4,sd=1)

標準偏差は、

sd(x1+x2)

## [1] 1.425891

平均は、

mean(x1+x2)

## [1] 5.995558

2.2 差の確率分布

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)

## No id variables; using all as measure variables

g <- ggplot(df, aes(x = value, fill = variable))
g <- g + geom_histogram(position = "identity", alpha = 0.8)
g <- g + scale_fill_npg()
plot(g)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

標準偏差は、

sd(x1-x2)

## [1] 1.411946

平均は、

mean(x1-x2)

## [1] -2.046432

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)\)
が成立します。（証明については、ここでは立ち入りません）

前出の例では、標準偏差 sd=1 は、分散 1 を意味し、 \(V(x_1)=1,~V(x_2)=1\) から
　　　\(V(x_1+x_2)=V(x_1)+V(x_2)=2\) 　　　\(\therefore~~sd(x_1+x_2)=\sqrt{2}\)
と考えると良いようです。

3 探究

正規乱数を用いる方法は、様々な確率変数の確率分布を考える上で役に立ちます。次の確率分布を正規乱数を用いて考えましょう。

正規分布の積
正規分布の商
正規分布の平方和

3.1 カイ二乗分布

正規分布の平方和はカイ二乗分布に従うことが知られています。
どのような確率分布になるのか正規乱数を用いて調べてみましょう。
rnorm(1000)により、平均０、標準偏差１の正規乱数を1000個生成します。

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

独立な5つの標準正規分布の平方和は、自由度５のカイ二乗分布に従うことが確かめられます。 dchisq(x,5) は自由度５のカイ二乗分布を表す関数です。

3.2 様々な確率分布

このほかにも特徴のある確率分布があります。

ポアソン分布
指数分布
ガンマ分布

これらの確率分布について調べてみましょう。

確率変数と確率分布