生物統計学 第3回　連続確率変数

連続確率変数

確率密度関数、累積密度関数

連続確率変数は、実数\(R\)上の区間\((a,b)\)で値を取る確率密度関数（probability density function: PDF）\(f(x)\)により、その変数が特定される。PDFは、非負であり、積分すると1となる。連続確率変数\(X\)が区間\((a,b)\)に値をもつ確率は、 \[ P[X\in (a,b)]=\int_a^bf(x)dx \] である。累積密度関数（CDF）は、 \[ F(x)=P(X\le x)=\int_{-\infty}^xf(t)dt \] である。連続確率変数\(X\)が区間\((a,b)\)に値をもつ確率は、CDFを用いて \[ P[X\in (a,b)]=F(b)-F(a) \] と表現できる。

期待値、モーメント

\(X\)の期待値は、 \[ E(X)=\int_R xf(x)dx \] と表され、任意の確率変数の関数\(g(X)\)の期待値は、 \[ E[g(X)]=\int_Rg(x)f(x)dx \] と表される。

連続確率変数の\(k\)次モーメントは、 \[ E(X^k)=\int_Rx^k f(x)dx \] と表され、\(k\)次の中心モーメントは、 \[ E[(X-E(X))^k]=\int_R (x-E(X))^k f(x)dx \] と表される。

なお、離散確率変数の場合と同様、1次モーメントは期待値であり、2次の中心モーメントが分散である。また、離散確率変数の場合と同様、2次モーメントと1次モーメントを用いて、分散を以下のように計算することもできる。

\[ Var(X) = E(X^2) - E(X)^2 \] なお、分散のもつ性質は、離散確率変数の場合と同じ。 すなわち、\(a, b\)を任意の実数、\(X\)を確率変数とするとき、 \[ Var(aX) = a^2Var(X) \\ Var(b) = 0 \\ Var(aX + b) = a^2Var(X) \]

積率母関数

　連続変数\(X\)の積率母関数（moment-generating function）は、 \[ m(t)=E(e^{tX})=\int_R e^{tx} f(x)dx \] である。ここで、 \[ m^{(k)} (t)=\int_R x^k e^{tx} f(x)dx \] であることから、 \[ E(X^k)=m^{(k)} (0) \] となり、離散確率変数の場合と同様、積率母関数の\(k\)次微分を行い、\(t=0\)とすることで\(k\)次モーメントを求めることができる。

幾つかの代表的な連続分布

一様分布

確率変数\(X\)は、以下のPDFをもつ場合に、一様分布に従う。

\[ f_X (x)=\Bigl\{ \begin{array}{ll} \frac{1}{b-a} & a\le x\le b \\ 0 & それ以外 \end{array} \] XのCDFは、 \[ F_X (x)=\Biggl\{ \begin{array}{ll} 0 & x < a \\ \frac{x-a}{b-a} & a \le x \le b \\ 1 & x > b \end{array} \] である。

\(X\)の積率母関数（付録1）は、 \[ m(t)=\frac{e^{tb}-e^{ta}}{t(b-a)} \] で、\(k\)次モーメントは、 \[ E(X^k)=\frac{1}{k+1} \sum_{i=0}^ka^i b^{k-i} \] である。

期待値と分散は、 \[ E(X)=\frac{a+b}{2} \] \[ Var(X)=E(X^2)-E(X)^2 = \frac{b^2+ab+a^2}{3} - \frac{b^2 + 2ab + a^2}{4} = \frac{(b-a)^2}{12} \] となる。

なお、\(a=0,b=1\)のとき\(X\)の分布を標準一様（standard uniform）分布とよぶ。また、2つの独立した標準一様分布の和は、三角形の密度分布 \[ f_X (x)=\Biggl\{ \begin{array}{ll} x & 0 \le x \le 1 \\ 2-x & 1 \le x \le 2 \\ 0 & それ以外 \end{array} \] となり、「魔女の帽子」ともよばれる分布となる。

op <- par(mfrow = c(1,2)) # 図を1行2列で表示する準備
# 魔女の帽子
x <- runif(10000)  # 一様分布から1万個の乱数を発生
hist(x)
x <- x + runif(10000)  # 独立にさらに1万個の乱数を発生
hist(x)

par(op) # 設定を基に戻す（重要！)

指数分布

指数確率変数の確率密度関数（PDF）は、 \[ f_X (x)=\Bigl\{ \begin{array}{ll} \lambda e^{-\lambda x} & x \ge 0 \\ 0 & それ以外 \end{array} \] ここで、\(\lambda>0\)はrateパラメータとよばれる。累積密度関数（CDF）は、 \[ F(x)=1-e^{-\lambda x} \] である。

\(X\)の積率母関数は、 \[ m(t)=\frac{\lambda}{\lambda-t} \] （ただし、\(t< \lambda\)）であり、\(k\)次モーメントは、 \[ E(X^k)=\frac{k!}{\lambda^k} \] である。

　期待値と分散は、それぞれ、 \[ E(X)=\frac{1}{\lambda} \] \[ Var(X)=E(X^2) - E(X)^2 = \frac{2}{\lambda^2} - \frac{1}{\lambda^2} = \frac{1}{\lambda^2} \] となる。

　指数分布は、ポワソン分布と重要な関連がある。ある事象が生じる時間間隔\(X\)が指数分布に従うとき、\(t\)時間内に生じる事象の回数\(n\)はポワソン分布に従う。より一般的には、\(X_1,X_2,\dots ,X_n\)がそれぞれ同じ指数分布に従うときに、それらの和\(S_n=X_1+X_2+\cdots +X_n\)を定義すると、全ての正の値\(t\)に対して、 \[ P(S_n\le t<S_{n+1})=p_Y(n) \] が成り立つ（左辺は\(t\)時間内に注目している事象がちょうど\(n\)回起こる確率）。ここで、\(p_Y (n)\)は、パラメータ\(\lambda t\)のポアソン確率変数\(Y\)の関数である。

　なお、指数分布は、無記憶性とよばれる性質をもつ。すなわち、 \[ P(X\le u+v │ X > u)=P(X\le v) \]

これは、 \[ P(X\le u+v │ X > u)=\frac{P(X\le u+v)-P(X\le u)}{1-P(X\le u)}\\ = \frac{F(u+v) - F(u)}{1-F(u)} = \frac{(1-e^{-\lambda(u+v)}) - (1-e^{-\lambda u})}{1-(1-e^{-\lambda u})} \\ = \frac{e^{-\lambda u} \cdot e^{-\lambda v} - e^{-\lambda u}}{e^{-\lambda u}} = 1-e^{-\lambda v} = P(X \le v) \] として導出できる。

例

ある農業機械は、購入後5年以内に30％が故障して使えなくなる。この機械が使えなくなるまでの時間\(T\)が、未知のパラメータ\(\lambda\)をもつ指数分布に従うとき、購入後10年以上使うことができる確率を求めよ。

　まず最初に、生存率から指数分布のパラメータ\(\lambda\)を求める。 \[ P(T>t)=1-F(t)=\exp(-\lambda t) \] であることから、\(P(T>5)=exp(-5\lambda)=1 - 0.3\)。したがって、 \[ \lambda=-\frac{1}{5}\log(0.7) \]

次に、10年間以上利用できる確率をCDFを用いて計算する。 \[ P(T>10)=1-F(10)=\exp⁡(-10\lambda) \]

具体的には、

lamda <- - 1/5 * log(0.7)
lamda

## [1] 0.07133499

exp(-10 * lamda)

## [1] 0.49

# Rの間数を用いて計算する
1 - pexp(10, lamda)

## [1] 0.49

グラフを描いてみる。

# draw graph
x <- 1:20
plot(x, y = 1 - pexp(x, lamda), ylim = c(0,1), type = "l")

正規分布（normal distribution）

正規分布のPDFは、以下のように表される。 \[ f(x)=\frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right) \]

正規分布のCDFは次のように表される。 \[ F(x)=\int_{-\infty}^x \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left(-\frac{(x-\mu)^2}{2\sigma}\right) \] となる。

正規分布の積率母関数（付録1）は、 \[ m_X(t) = \exp\left(\mu t + \frac{\sigma^2t^2}{2}\right) \] である。

期待値と分散（付録1）は、 \[ E(X)=\mu \] \[ Var(X)=\sigma^2 \] となる。

標準正規分布と確率計算

　正規分布は、その平均と標準偏差が変わると分布の位置や形状が変化する。言い換えると、平均と標準偏差が決まると分布の形も決まる。確率変数\(X\)が、平均\(\mu\)、標準偏差\(\sigma\)（あるいは分散\(\sigma^2\)）の正規分布に従っていることを、以下のように書き表す。 \[ X \sim N(\mu,\sigma^2) \] 例えば、平均50, 標準偏差10の正規分布は、\(N(50,10^2)\)あるいは\(N(50,100)\)と表される。

　正規分布では平均と標準偏差が変わると分布の位置や形状が変化するが、確率変数に次のような変換（標準化）を施すと、変換後の変数は平均0、分散1の正規分布（この正規分布を、標準正規分布（standard normal distribution）とよぶ）に従う。すなわち、 \[ X\sim N(\mu,\sigma^2) \longrightarrow Z=(X-\mu)/\sigma \sim N(0,1) \]

　こうして変換しておけば、例えば、積分計算が必要になるCDFを標準正規分布について予め計算しておき、それを他の正規分布に利用することができる。なお、標準正規分布のPDFは小文字の\(\phi\)を使って \[ \phi(x)=\frac{1}{\sqrt{2\pi}} \exp\left(-\frac{x^2}{2}\right)⁡ \] と表され、CDFは大文字の\(\Phi\)を使って \[ \Phi(x)=\int_{-\infty}^x \phi(x)dx=\int_{-\infty}^x \frac{1}{\sqrt{2\pi}}\exp\left(-\frac{x^2}{2}\right) \] と表される。

例

ある販売用の果物の重さが平均80g、標準偏差20ｇ正規分布\(N(80,20^2)\)に従うとする。このとき、

1. 40g以下の果実は売り物にならないとすると、全体の何%が販売可能な果実か
2. 販売できる果実のうち60g以下の果実を安売りにまわすとすると、全体の何%が安売り品になるか
3. 60g以上の品については、重い方からA, B, Cの3等級に分け、それぞれの割合を25, 50, 25%とするには、等級AとB、等級BとCの境界をどのように決めればよいか

　まず、果実の重さを確率変数\(X\)として、先ほどの変換を用いると、以下のようになる。すなわち、 \[ P(40<X)=P\left(\frac{40-80}{20}<Z\right)=1-P\left(Z\le \frac{40-80}{20}\right)=1-\Phi\left(\frac{40-80}{20}\right) \] これは、Rを用いて簡単に計算することができる。

1 - pnorm((40 - 80) / 20)

## [1] 0.9772499

計算の結果、果実の97.7%を販売できることが分かる。

　次に、安売りにまわされる果実の割合も同様に計算できる。すなわち、 \[ P(40<X≤60)=P\left(\frac{40-80}{20}<Z\le \frac{60-80}{20}\right)=\Phi\left(\frac{60-80}{20}\right) - \Phi\left(\frac{40-80}{20}\right) \] これを、Rを用いて計算すると、

pnorm((60 - 80) / 20) - pnorm((40 - 80) / 20)

## [1] 0.1359051

結果、13.6%の果実が安売りにまわることが分かる。

　最後に、等級分けをするための境界を求めてみる。まず、等級AとBの境界を求める。等級Aに属するのは、果実重が60gよりも重い果実の25%となる。従って、その割合\(p_A\)は、 \[ p_A=\frac{1}{4} P(60<X)=\frac{1}{4} \left[1-\Phi\left(\frac{60-80}{20}\right)\right] \] となる。いま、等級AとBの境界の重さを\(a\)とおくと、 \[ P(a<X)=1-\Phi \left(\frac{a-80}{20}\right)=p_A \\ \leftrightarrow \Phi \left(\frac{a-80}{20}\right) = 1-p_A \] となる必要があることが分かる。Rでは、正規分布の累積分布関数の逆関数（累積確率から確率変数の値を求める関数）としてqnormがある。これを用いて上述した計算を行う。

# p_Aの計算と表示
pA <- 0.25 * (1 - pnorm((60 - 80) / 20))
pA

## [1] 0.2103362

# Zにおける境界の計算
zA <- qnorm(1 - pA)
zA

## [1] 0.8052553

# Xへの変換
xA <- 20 * zA + 80
xA

## [1] 96.10511

その結果、等級AとBの境界は96.1gとなる。

　同様に、等級BとCの境界を\(b\)とおくと、\(b\)は以下を満たす必要がある。 \[ \Phi\left(\frac{b-80}{20}\right)=1-\frac{3}{4} P(60<X)=1-3 p_A \] そこで、先に計算されている\(p_A\)も用いて計算すると、

# ZにおけるB-C境界の計算
zB <- qnorm(1 - 3 * pA)
zB

## [1] -0.3345257

# Xへの変換
xB <- 20 * zB + 80
xB

## [1] 73.30949

計算の結果、等級BとCの境界は73.3gとなることが分かる。

　最後に、\(N(80,20^2)\)に従う乱数を発生させて、10,000個の果実の重さを擬似的に生成して、その分類（廃棄、安売り、ABC等級付け）を行い、上述した計算結果と比較してみる。正規分布に従う乱数は、関数rnormを用いて発生することができる。

fw <- rnorm(10000, 80, 20)
# 販売できる果実の数と割合
sum(fw > 40)

## [1] 9781

sum(fw > 40) / length(fw)

## [1] 0.9781

# 安売りにまわされる果実の割合
sum(fw > 40 & fw <= 60) / length(fw)

## [1] 0.1341

# 安売りせずにすむ果実のうち等級A、B、Cの割合
sum(fw > xA) / sum(fw > 60)

## [1] 0.2441943

sum(fw > xB & fw <= xA) / sum(fw > 60)

## [1] 0.5014218

sum(fw > 60 & fw <= xB) / sum(fw > 60)

## [1] 0.2543839

以上をヒストグラムを描いて確認する。

# ヒストグラムを用いた図示
hist(fw)
abline(v = c(40, 60, xB, xA), col = c("red", "yellow", "green", "blue"))

正規分布の再生性

　正規分布の重要な特性のひとつに、再生性（reproductive property）という特性がある。これは、独立した正規分布に従う確率変数のあらゆる線形結合が、やはり、正規分布に従うというものである。正規分布では平均と分散が分布の形を決定しているため、線形結合に含まれる各正規分布の平均と分散さえわかれば、線形結合から得られる正規分布の形も規定できる。 　今、\(X_1,X_2,…,X_n\)を独立した正規分布に従う確率変数とし、\(X_i~N(\mu_i,\sigma_i^2)\)とする。このとき、いかなる線形結合の重み\(a_1,a_2,\dots,a_n\)に対しても以下が成り立つ。 \[ a_1 X_1+a_2 X_2+\cdots+a_n X_n = \sum_{i=1}^n a_i X_i \sim N(\mu,\sigma^2) \] ここで、 \[ \mu=a_1 \mu_1+a_2 \mu_2+\cdots+a_n \mu_n=\sum_{i=1}^n a_i \mu_i \\ \sigma^2=a_1^2 \sigma_1^2+a_2^2 \sigma_2^2+\cdots+a_n^2 \sigma_n^2=\sum_{i=1}^n a_i^2 \sigma_i^2 \] 　重要な特殊な場合として、次の2つの場合、(1) \(a_1=1,a_2=-1\)、(2) \(a_1=\cdots=a_n=1/n\)が挙げられる。(1)の場合は、2つの正規分布の差の分布を表しており、平均は、2つの正規分布の平均の差となるが、分散は、2つの正規分布の分散の「和」となることに注意する。 (2)の場合は、複数の正規分布の算術平均\(\bar X\) の分布となる。例えば、\(X_1,X_2,\dots,X_n\)がそれぞれ独立に同じ正規分布\(N(\mu,\sigma^2)\)に従っている場合は、その標本平均\(\bar X =(X_1+\cdots+X_n)/n\)は正規分布\(N(\mu,\sigma^2⁄n)\)に従う。したがって、\(X_i\)の分散と\(\bar X\) の分散には、以下の関係がある。 \[ \sigma_{\bar X}^2=\frac{\sigma^2}{n} \]

例

　それぞれが30aの広さをもつ2枚の田圃でイネの品種AとBを栽培することを考える。品種Aは、多収性は高いが少し不安定で、1枚の田圃あたりの収量は、平均1.7t、標準偏差0.3tの正規分布\(N(1.7, 0.3^2 )=N(1.7,0.09)\)に従う。いっぽう、品種Bは、多収性は品種Aに劣るが安定しており、1枚の田圃あたりの収量は、平均1.5t、標準偏差0.1tの正規分布\(N(1.5,0.1^2 )=N(1.5,0.01)\)に従う。また、2枚の田圃での収量は独立に各正規分布に従っていると仮定する。このとき、品種Aの収量が品種Bの収量を下回る確率はいくらになるか。また、品種AとBの標準偏差が逆（品種Aの標準偏差が0.1t、品種Bの標準偏差が0.3t）の場合には、結果はどのように変わるか。

　品種Aの収量\(X_A\)と品種Bの収量\(X_B\)を比較するには、収量の差\(D=X_A-X_B\)の分布を考えればよく、品種Aの収量が品種Bの収量を下回るのは、\(D\)が0未満になる場合になる。\(D\)が従う分布は、\(N(0.2,0.1)\)であり、\(D\)が0未満になる確率は、 \[ P(D<0)=\Phi\left(\frac{0-0.2}{\sqrt{0.1}}\right) \] これをRで計算すると

# 問題2(1)、品種Aが品種Bを下回る確率
pnorm((0 - 0.2) / sqrt(0.1))

## [1] 0.2635446

となり、26.4%の確率で品種Aが品種Bを下回ることが分かる。つまり、品種Aは平均すると品種Bよりも0.2t多く収穫できるものの、およそ4回に1回は品種Bの収量を下回る。なお、品種AとBの収量の分散を逆（品種Aの分散を\(0.1^2\)、品種Bの分散を\(0.3^2\)）にしても、\(D\)の従う正規分布の分散は変わらないため同じ結果になる。 　では、実際の値の分布を、100,000回のシミュレーションを行って眺めてみる。

n.sim <- 100000# 100000回のシミュレーション
# 品種Aの収量
a <- rnorm(n.sim, 1.7, 0.3)
# 品種Bの収量
b <- rnorm(n.sim, 1.5, 0.1)
# 品種AとBの収量差
d <- a - b
# 差が負になる確率
sum(d < 0) / length(d)

## [1] 0.26343

# ヒストグラムの描画
# 最初に値の範囲と区切り位置だけを決める
tmp <- hist(c(a, b, d), plot = F)
# 順次描画
par(mfrow = c(3, 1))
hist(a, breaks = tmp$breaks, main = "Var A")
hist(b, breaks = tmp$breaks, main = "Var B")
hist(d, breaks = tmp$breaks, main = "A - B")

# 描画設定を元に戻す
par(mfrow = c(1, 1))

中心極限定理

　中心極限定理（central limit theorem: CLM）は、正規分布を、様々な確率分布の中でも、最も重要で根幹となるものとする定理である。前節で、独立した正規分布の線形結合は、また、正規分布となるという再生性について説明した。すなわち、\(X_1,X_2,\dots,X_n\)がそれぞれ独立に\(N(\mu,\sigma^2)\)に従っている場合、その和の分布は、 \[ \sum_{i=1}^n X_i \sim N(n\mu,n\sigma^2) \] また、平均の分布は、 \[ \bar X = \frac{1}{n} \sum_{i=1}^n X_i \sim N\left(\mu,\frac{\sigma^2}{n}\right) \] と表される。

　中心極限定理は、\(\sum_{i=1}^n X_i\)あるいは\(\bar X\) が「近似的に」正規分布にあてはまるためには、\(X_1,X_2,\dots,X_n\)が必ずしも正規分布に従っている必要はない、という定理である。正規分布への近似の程度は、\(n\)の大きさが30を超えていれば非常に良いことが知られている。なお、中心極限定理が成り立つためには、\(X_1,X_2,\dots,X_n\)は正規分布に従っている必要はないが、\(X_i\)が互いに独立で、同じ分布に従い、かつ、有限の分散と平均をもつ、という条件が満たされているい必要がある。逆に、この条件が満たされていれば、どのような場合でも中心極限定理が近似的に成り立つ。

中心極限定理： 　\(X_1,X_2,\dots,X_n\)がそれぞれ独立に、同じ分布に従う確率変数で、その分布の平均と分散が有限の場合、 \[ \sum_{i=1}^n X_i \sim N(n\mu,n\sigma^2) \] および \[ \bar X = \frac{1}{n} \sum_{i=1}^n X_i \sim N\left(\mu,\frac{\sigma^2}{n}\right) \] が近似的に成り立つ。

では、中心極限定理を、Rを用いたシミュレーションで確認してみる。ここでは、\(X_1,X_2,\dots,X_n\)がそれぞれ独立に、\([0,1]\)の範囲の一様分布に従うと仮定する。\(n\)を1〜6まで変化させ、\(\sum_{i=1}^n X_i\)を計算してヒストグラムを作成してみる。なお、\([a,b]\)の範囲の一様分布は、平均と分散が \[ \frac{b – a}{2} \\ \frac{(b – a)^2}{12} \] に従う。そこから中心極限定理から期待される正規分布の平均と分散を計算して、重ね描きしてみる。

# 中心極限定理のシミュレーション
# 一様乱数からの標本を600000個準備して、100000×6の行列に収める
n.sim <- 100000 # 100000回のシミュレーション
n <- 6 # 足し合わせる変数の最大数
r <- runif(n.sim * n) # 600000個の一様乱数
d <- matrix(r, nrow = n.sim) # 100000行（6列)の行列に乱数を収める
# 予め6個分の図のスペースを準備
par(mfrow = c(2, 3))
# 1列目（1変数目）だけを用いてヒストグラムを描画
x1 <- d[, 1] 
hist(x1, freq = F)
# 中心極限定理から期待される分布曲線の描画
x <- seq(0, 1, 0.01)
y <- dnorm(x, 1/2, sqrt(1/12))
lines(x, y, col = "green")
# 準備しておいた行列から1列目と2列目の和を計算
x1_2 <- d[, 1] + d[, 2]
# ヒストグラムの描画
hist(x1_2, freq = F)
# 中心極限定理から期待される分布曲線の描画
x <- seq(0, 2, 0.01)
y <- dnorm(x, 2/2, sqrt(2/12))
lines(x, y, col = "green")
# 以下、3...6まで繰り返す。
# 行和をとるのにapply変数を使うと便利である
x1_3 <- apply(d[, 1:3], 1, sum) # 行列の1-3列の和をとる
hist(x1_3, freq = F)
x <- seq(0, 3, 0.01)
y <- dnorm(x, 3/2, sqrt(3/12))
lines(x, y, col = "green")
x1_4 <- apply(d[, 1:4], 1, sum) # 行列の1-4列の和をとる
hist(x1_4, freq = F)
x <- seq(0, 4, 0.01)
y <- dnorm(x, 4/2, sqrt(4/12))
lines(x, y, col = "green")
x1_5 <- apply(d[, 1:5], 1, sum) # 行列の1-5列の和をとる
hist(x1_5, freq = F)
x <- seq(0, 5, 0.01)
y <- dnorm(x, 5/2, sqrt(5/12))
lines(x, y, col = "green")
x1_6 <- apply(d[, 1:6], 1, sum) # 行列の1-6列の和をとる
hist(x1_6, freq = F)
x <- seq(0, 6, 0.01)
y <- dnorm(x, 6/2, sqrt(6/12))
lines(x, y, col = "green")

# 図のスペースを元に戻す
par(mfrow = c(1, 1))

正規分布に関連する分布

　ここで、正規分布に関連する4つの重要な分布、カイ二乗分布、スチューデントのt分布、F分布、対数正規分布、について、簡単に紹介する。これら分布と正規分布の関係は、独立した標準正規分布の関数のかたちで表される。ここで、\(Z_1,Z_2,\dots,Z_n\)を\(n\)個の独立した標準正規分布に従う確率変数であるとする。このとき、4つの確率分布と正規分布の関係は以下のように表せる。

カイ二乗（Chi-square）分布 - \(Z_i\)の平方和、すなわち、\(Z_1^2+\cdots+Z_n^2\)は、自由度\(n\)のカイ二乗分布\(\chi_n^2\)に従う。すなわち、 \[ \chi_n^2 \sim Z_1^2+\cdots+Z_n^2 \]

スチューデントの\(t\)（Student’s \(t\)）分布 - 標準正規分布に従う確率変数\(Z\)と、それとは独立したカイ二乗分布に従う確率変数\(\chi^2\)をその自由度で除算したものの平方根との比は、自由度\(n\)のスチューデントの\(t\)分布\(t_n\)に従う。すなわち、 \[ t_n\sim \frac{Z}{\sqrt{\frac{\chi_n^2}{n}}} \]

\(F\)分布 - 2つの独立したカイ二乗分布に従う確率変数を、それぞれの自由度\(m, n\)で除算したものの比は、自由度\(m, n\)の\(F\)分布\(F_{m,n}\)に従う。すなわち、 \[ F_{m,n}\sim \frac{\chi_m^2⁄m}{\chi_n^2⁄n} \]

対数正規（log-normal）分布 - 対数正規分布とは、その確率変数の対数が正規分布に従っているものをいう。すなわち、\(X\)が対数正規分布に従っている場合、\(Y=\log X\)が正規分布に従う。

　カイ二乗分布と正規分布の関係をシミュレーションで確認してみる。ここでは、標準正規分布にしたがう10個の独立な確率変数を発生し、この2乗和を計算し、自由度10のカイ二乗分布と分布の形を比較してみる。 　実際の計算は、100000×10個の正規分布を発生させ、その2乗を100000行×10列の行列として整形し、各行について和をとるという作業を行って、上述した試行を100000回行う。

# カイ二乗分布と正規分布の関係を調べるシミュレーション
n.sim <- 10000
n <- 10
r <- matrix(rnorm(n.sim * n)^2, nrow = n.sim)
a <- apply(r, 1, sum)
hist(a, freq = F)
x <- seq(0, 40, 0.1)
y <- dchisq(x, n)
lines(x, y, col = "green")

　次に、t分布に従う乱数を発生させ、そのヒストグラムと正規分布の形状を比較してみる。関数\(rt\)を用いて自由度3の\(t\)分布に従う乱数を10,000個発生させている。そのヒストグラムを描き、その上に、発生させた乱数と同じ平均、分散をもつ正規分布と、自由度3の\(t\)分布を重ね描きしてみる。

# t分布に従う乱数の分布と正規分布の比較
# 自由度3のt分布からの乱数を100000個発生
n.sim <- 100000
d <- rt(n.sim, 3)
# ヒストグラムの描画。1ずつでヒストグラムの区間を区切っている
hist(d, freq = F, breaks = seq(floor(min(d)), ceiling(max(d)), 1))
# 同じ平均と標準偏差をもつ正規分布を描画
x <- seq(floor(min(d)), ceiling(max(d)), 0.1)
y <- dnorm(x, mean(d), sd(d))
lines(x, y, col = "red")
# 自由度3のt分布を描画
y <- dt(x, 3)
lines(x, y, col = "green")

　最後に、\(F\)分布をシミュレーションで作成して、ヒストグラムを描いてみる。

n.sim <- 100000
na <- 8
r <- matrix(rnorm(n.sim * na)^2, nrow = n.sim)
a <- apply(r, 1, sum)
nb <- 10
r <- matrix(rnorm(n.sim * nb)^2, nrow = n.sim)
b <- apply(r, 1, sum)
d <- (a / na) / (b / nb)
hist(d, freq = F, breaks = seq(0, ceiling(max(d)), 0.1))
x <- seq(0, max(d), 0.01)
y <- df(x, na, nb)
lines(x, y, col = "green")

付録

付録1 正規分布の積率母関数、期待値、分散

正規分布の積率母関数 \[ m_X(t) = E(e^{tX}) \\ = \int_{-\infty}^\infty e^{tX} f(x) dx \\ = \int_{-\infty}^\infty \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left[-\frac{(x-\mu)^2}{2\sigma^2} + tx\right] dx \\ = \int_{-\infty}^\infty \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left[-\frac{1}{2\sigma^2}(x^2 -2\mu x + \mu^2 - 2\sigma^2 tx)\right] dx \\ = \int_{-\infty}^\infty \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left[-\frac{1}{2\sigma^2}\left((x -(\mu + \sigma^2 t))^2 + 2\mu \sigma^2 t + \sigma^4 t^2\right)\right] dx \\ = \int_{-\infty}^\infty \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left[-\frac{(x -(\mu + \sigma^2 t))^2}{2\sigma^2}+ \mu t + \frac{\sigma^2 t^2}{2}\right] dx \\ = \exp{\left(\mu t + \frac{\sigma^2 t^2}{2}\right)}\int_{-\infty}^\infty \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left[-\frac{(x -(\mu + \sigma^2 t))^2}{2\sigma^2}\right] dx \] ここで、 \[ \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left[-\frac{(x -(\mu + \sigma^2 t))^2}{2\sigma^2}\right] \] は、期待値\(\mu + \sigma^2 t\)、分散\(\sigma^2\)の正規分布の密度関数と考えられるので、 \[ \int_{-\infty}^\infty \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left[-\frac{(x -(\mu + \sigma^2 t))^2}{2\sigma^2}\right] dx = 1 \] したがって、 \[ m_X(t) = \exp{\left(\mu t + \frac{\sigma^2 t^2}{2}\right)} \]

1次モーメントは、 \[ m'_X(t) = \left(\mu t + \frac{\sigma^2 t^2}{2}\right)' \exp{\left(\mu t + \frac{\sigma^2 t^2}{2}\right)} \\ = (\mu + \sigma^2 t) \exp{\left(\mu t + \frac{\sigma^2 t^2}{2}\right)}\\ \therefore m'_X(0) = \mu \]

したがって、期待値は、 \[ E(X) = m'_X(0) = \mu \] となる。

2次モーメントは、 \[ m''_X(t) = \left[(\mu + \sigma^2 t) \exp{\left(\mu t + \frac{\sigma^2 t^2}{2}\right)}\right]' \\ = (\mu + \sigma^2 t)' \exp{\left(\mu t + \frac{\sigma^2 t^2}{2}\right)} + (\mu + \sigma^2 t)^2 \exp{\left(\mu t + \frac{\sigma^2 t^2}{2}\right)}\\ = \left(\sigma^2 + (\mu + \sigma^2 t)^2\right)\exp{\left(\mu t + \frac{\sigma^2 t^2}{2}\right)} \\ \therefore m''_X(0) = \sigma^2 + \mu^2 \]

したがって、分散は、 \[ var(X) = E(X^2) - E(X)^2 = (\sigma^2 + \mu^2) - \mu^2 = \sigma^2 \]