演習問題

標準偏差 … データのばらつきを議論したいときに利用
標準誤差 … 真の平均の不偏推定値 (\(\mu = \frac{1}{N}\sum_i x_i\)) のばらつきを議論したいときに利用

前の講義資料を参照してください。

統計の概念図 (再掲)

需要な点なので、再掲する。

母集団からランダムにサンプリングされるデータを確率変数とみなして、計測されるデータから未知の母集団の平均や分散を推定する。これが統計の枠組みである。
具体的に計測されたデータ (e.g., 150cm, 165cm, … などの身長データ) の平均を計算すると、期待値としては真の平均に一致することを示した。ただし、あくまで期待値として一致するのみであり、分散 \(\frac{\sigma^2}{N}\) と伴いばらつく (\(\sigma^2\)も未知)。

大数の法則 (再び)

データ数 \(N\)、データ \(x_1, ..., x_N\) のとき、データ平均 \(\bar{x}\) は、 \[\begin{equation} \bar{x} = \frac{1}{N} \sum_{i = 1}^N x_i \end{equation}\] となる。このデータ平均の母集団分布での期待値 \(E[\cdot]\) は、未知の母集団平均 \(\mu\) と用いて \[\begin{equation} E[\bar{x}] = \frac{1}{N} \sum_{i = 1}^N E[x_i] = \mu \end{equation}\] となる。つまり、データ平均は未知の母集団平均 \(\mu\) の不偏推定量である。

加えて、データ平均の分散の期待値 \(E[\left( \frac{1}{N}\sum_{i=1}^Nx_i - \mu \right)^2]\) は、各データが独立であることを利用すると、未知の母集団分散 \(\sigma^2\) と用いて \[\begin{align} \left( \frac{1}{N}\sum_{i=1}^Nx_i - \mu \right)^2 &= \left( \frac{1}{N}\sum_{i=1}^N (x_i - \mu) \right)^2 \nonumber\\ &= \frac{1}{N^2} \sum_{i=1}^N E[(x_i - \mu)^2] \nonumber\\ &= \frac{1}{N} \sigma^2 \end{align}\] となる。すなわち、 \(N \to \infty\) において、データ平均のばらつきは 0 に収束していき、\(\mu\) に一致する。

前回の講義資料にもあったように、平均や分散が変数変換により受ける影響は \[\begin{equation} E[aX + b] = aE[X] + b \end{equation}\] \[\begin{equation} \mathrm{Var}[aX + b] = a^2\mathrm{Var}[X] \end{equation}\] となっており、これを利用するとサクッと計算できる。

以上より、\(N \to \infty\) において、データ平均を計算すれば確率1で\(\mu\) に一致することが示せた。このような法則を、大数の法則と呼ぶ。

この大数の法則を数値的に確かめてみよう。下記は、N = 1 (赤)、N = 5 (緑)、N = 100 (青)、N = 10000 (黒) にてデータ平均を計算したシミュレーション結果である。このデータ平均の計算を1000回繰り返している。真の平均は0である(本来は未知のものだがシミュレーションなのでわかる)。

N1 = 1 # 10サンプル
N2 = 5
N3 = 10
N4 = 100
N5 = 10000
Ite = 1000 # 1000回繰り返す
  
mu1 = matrix(0, 1, Ite)
mu2 = matrix(0, 1, Ite)
mu3 = matrix(0, 1, Ite)
mu4 = matrix(0, 1, Ite)
mu5 = matrix(0, 1, Ite)

for(i in 1:Ite){
  dat = rnorm(N5, 0, 1)
  mu1[i] = mean(dat[1:N1])

  mu2[i] = mean(dat[1:N2])

  mu3[i] = mean(dat[1:N3])

  mu4[i] = mean(dat[1:N4])
  
  mu5[i] = mean(dat)
}

hist(mu1, prob = T, ylim = c(0, 1), col="#ff000020")
hist(mu2, prob = T, ylim = c(0, 1), col="#00ff0020", add = T)
# hist(mu3, prob = T, ylim = c(0, 1), col="#0000ff20", add = T)
hist(mu4, prob = T, ylim = c(0, 1), col="#0000ff20", add = T)
hist(mu5, prob = T, ylim = c(0, 1), col="#00000020", add = T)

上の図を見て分かる通り、Nが大きくなるにつれて、標本分布が真の平均 0 に近づいていく様子が見てとれると思う。 分布が近づくとは、ここでは平均が0かつ分散が小さくなることを意味する。

以上より、大数の法則を数値的に確かめることができた。

標本分散 (データ分散)

データ平均が真の母集団平均の不偏推定量であるからに、データ分散も真の母集団分散の不偏推定量ではないかと考えられる。確かめてみよう。

N1 = 2 # N=1では分散は常に0になる
N2 = 5
N3 = 10
N4 = 100
N5 = 10000
Ite = 1000 # 1000回繰り返す
  
mu1 = matrix(0, 1, Ite)
mu2 = matrix(0, 1, Ite)
mu3 = matrix(0, 1, Ite)
mu4 = matrix(0, 1, Ite)
mu5 = matrix(0, 1, Ite)

for(i in 1:Ite){
  dat = rnorm(N5, 0, 1)
  mu1[i] = mean(dat[1:N1]^2) - mean(dat[1:N1])^2

  mu2[i] = mean(dat[1:N2]^2) - mean(dat[1:N2])^2

  mu3[i] = mean(dat[1:N3]^2) - mean(dat[1:N3])^2

  mu4[i] = mean(dat[1:N4]^2) - mean(dat[1:N4])^2
  
  mu5[i] = mean(dat^2) - mean(dat)^2
}

t1 = c(N1, N2, N3, N4)
t2 = c(mean(mu1), mean(mu2), mean(mu3), mean(mu4))

print(t1)
## [1]   2   5  10 100
print(t2)
## [1] 0.5350531 0.8387085 0.9233865 0.9943879

真の母集団分散は1 (シミュレーションなのでわかる) に対して、N = 100のときはおおよそ1に一致している。 しかしながら、N = 2, 5, 10のときは、明らかに1より小さいことがわかる。

上のシミュレーションでは、2乗の平均 - 平均の2乗で分散を計算したが、少しこれを補正してみよう。

N1 = 2 # N=1では分散は常に0になる
N2 = 5
N3 = 10
N4 = 100
N5 = 10000
Ite = 1000 # 1000回繰り返す
  
mu1 = matrix(0, 1, Ite)
mu2 = matrix(0, 1, Ite)
mu3 = matrix(0, 1, Ite)
mu4 = matrix(0, 1, Ite)
mu5 = matrix(0, 1, Ite)

for(i in 1:Ite){
  dat = rnorm(N5, 0, 1)
  mu1[i] = mean(dat[1:N1]^2) - mean(dat[1:N1])^2

  mu2[i] = mean(dat[1:N2]^2) - mean(dat[1:N2])^2

  mu3[i] = mean(dat[1:N3]^2) - mean(dat[1:N3])^2

  mu4[i] = mean(dat[1:N4]^2) - mean(dat[1:N4])^2
  
  mu5[i] = mean(dat^2) - mean(dat)^2
}

# 補正する
mu1 = mu1*N1/(N1-1)
mu2 = mu2*N2/(N2-1)
mu3 = mu3*N3/(N3-1)
mu4 = mu4*N4/(N4-1)

t1 = c(N1, N2, N3, N4)
t2 = c(mean(mu1), mean(mu2), mean(mu3), mean(mu4))

print(t1)
## [1]   2   5  10 100
print(t2)
## [1] 1.0239324 1.0197106 1.0104988 0.9992491

今度は1に近くなっていることがわかる。ここでは、観測データ数を \(N\) として、 \[\begin{align} \bar{\mathrm{Var}}(x) &= \frac{N}{N-1} \frac{1}{N}\sum_i \left( x_i - \frac{1}{N} \sum_i x_i \right)^2 \nonumber\\ &= \frac{1}{N-1}\sum_i \left( x_i - \frac{1}{N} \sum_i x_i \right)^2 \end{align}\] として計算している。謎の係数 \(\frac{N}{N-1}\) が掛け算されていることに注意。

ちなみに、Rで分散を計算する標準コマンド var を利用すると、

N1 = 2 # N=1では分散は常に0になる
N2 = 5
N3 = 10
N4 = 100
N5 = 10000
Ite = 1000 # 1000回繰り返す
  
mu1 = matrix(0, 1, Ite)
mu2 = matrix(0, 1, Ite)
mu3 = matrix(0, 1, Ite)
mu4 = matrix(0, 1, Ite)
mu5 = matrix(0, 1, Ite)

for(i in 1:Ite){
  dat = rnorm(N5, 0, 1)
  mu1[i] = var(dat[1:N1])

  mu2[i] = var(dat[1:N2])

  mu3[i] = var(dat[1:N3])

  mu4[i] = var(dat[1:N4])
  
  mu5[i] = var(dat)
}

t1 = c(N1, N2, N3, N4)
t2 = c(mean(mu1), mean(mu2), mean(mu3), mean(mu4))

print(t1)
## [1]   2   5  10 100
print(t2)
## [1] 0.9840536 0.9846953 0.9867457 0.9992569

おおよそ 1 に近いことがわかる。つまり、Rでは標準で補正していることを意味する。

なぜ補正が必要なのか、計算してみよう (少しややこしい)。真の母集団平均を \(\mu\) (未知) 、母集団分散を \(\sigma^2\) (未知) として、 データ分散の期待値を計算する。

\[\begin{align} E\left[ \frac{1}{N}\sum_i \left( x_i - \frac{1}{N} \sum_i x_i \right)^2 \right] &= E\left[ \frac{1}{N}\sum_i \left( x_i - \mu + \mu - \frac{1}{N} \sum_i x_i \right)^2 \right] \nonumber\\ &= E\left[ \frac{1}{N}\sum_i \left( x_i - \mu \right)^2 + \frac{1}{N}\sum_i \left(\mu - \frac{1}{N} \sum_j x_j\right)^2 + 2\frac{1}{N}\sum_i\left( x_i - \mu \right)\left(\mu - \frac{1}{N} \sum_j x_j\right) \right] \end{align}\]

ここで、\(E\left[ \frac{1}{N}\sum_i \left( x_i - \mu \right)^2 \right] = \sigma^2\) (定義より)、\(E\left[\frac{1}{N}\sum_i \left(\mu - \frac{1}{N} \sum_j x_j\right)^2\right] = \frac{1}{N}\sum_i \frac{\sigma^2}{N}\) (データ平均と母集団平均の関係性から)となり、計算するべき値は残るは \(E\left[2\frac{1}{N}\sum_i\left( x_i - \mu \right)\left(\mu - \frac{1}{N} \sum_j x_j\right)\right]\)となる。丁寧に計算していく。

\[\begin{align} E\left[2\frac{1}{N}\sum_i\left( x_i - \mu \right)\left(\mu - \frac{1}{N} \sum_j x_j\right)\right] &= E\left[ 2\frac{1}{N}\sum_i \left( x_i\mu - \mu^2 - x_i \frac{1}{N} \sum_j x_j + \mu \frac{1}{N} \sum_j x_j\right) \right] \nonumber\\ &= E\left[ 2\frac{1}{N}\sum_i \left( - \frac{1}{N}x_i^2 - x_i \frac{1}{N} \sum_{j\neq i} x_j\right) \right] + 2\mu^2 \nonumber\\ &= -\frac{2}{N}(\mu^2 + \sigma^2) -2 \frac{N-1}{N}\mu^2 + 2\mu^2 \nonumber\\ &= -\frac{2}{N}\sigma^2 \end{align}\]

以上より、 \[\begin{align} E\left[ \frac{1}{N}\sum_i \left( x_i - \frac{1}{N} \sum_i x_i \right)^2 \right] &= \sigma^2 + \frac{\sigma^2}{N} -\frac{2}{N}\sigma^2 = \frac{N - 1}{N}\sigma^2 \end{align}\] を得る。すなわち、データ分散は未知の母集団分散の不偏推定量ではないことを意味する。
同様に、未知の母集団分散の不偏推定量は、 \[\begin{equation} \bar{\sigma}^2 = \frac{1}{N-1}\sum_i \left( x_i - \frac{1}{N} \sum_i x_i \right)^2 \end{equation}\] であり、確かに \(E[\bar{\sigma}^2] = \sigma^2\) を満たす。

少し不思議な気もするが、真の母集団平均の不偏推定量はデータ平均、真の母集団分散の不偏推定量は \(N-1\) で割り算をしたデータ分散となる。通常多くのソフトウェアでは、分散を計算するコマンドは \(N-1\) で割ることを実装していることが多いように見受けられる (そのため、分散の計算式を手打ちすると、ソフトウェアが出力した分散と異なるため困惑することがある)。

標準化

平均や分散を推定する場合には、不偏推定量を利用すれば大抵の場合には事足りる。しかしながら、これまた多くの場合には、平均や分散のみならずに回帰分析 \(y = f(a_1x_1 + a_2x_2 + ... + b)\) を行い、データ \(x_1, ..., x_N\) と \(y\) の関係性を検証することも多々ある (この講義では触れないが、回帰分析を含む一般化線形モデルは、分散分析と呼ばれる最も利用される統計手法の一つの拡張版である。回帰分析は諸刃の剣だが、メリット・デメリット理解して使いこなすと強力な武器になる)。

回帰分析など、機械学習や信号処理において、解析する前のデータ前処理は非常に重要である。その中でも重要なものの一つに、標準化がある。標準化とは、平均 \(\mu\) 、分散 \(\sigma^2\) の確率変数 \(X\) を、 \[\begin{equation} Z = \frac{X - \mu}{\sigma} \end{equation}\] として変換することである。この変換の意味を見ていく。

\(Z\)の期待値を計算してみると、 \[\begin{equation} E[Z] = \frac{E[X] - \mu}{\sigma} = 0 \end{equation}\] となり、分散は \[\begin{equation} \mathrm{Var}[Z] = \frac{1}{\sigma^2}\mathrm{Var}[X] = 1 \end{equation}\] となる。

すなわち、標準化により、確率変数 \(X\) を平均 0 、標準偏差 1 に変換することができる。 この標準化はしばしば、異なるデータ \(x_1\), \(x_2\), … を公平に比較する際や同時に解析する際に利用する。

中心極限定理 (かなり重要)

\(X_1, ..., X_n\)を確率分布から得られる、独立な確率変数とする。そして、各々平均 \(\mu\) 、分散 \(\sigma^2\) をもつものとする。このとき、
\[\begin{equation} S_n = X_1 + X_2 + ... + X_n \end{equation}\] とする。この \(S_n\)の平均は \(n\mu\)、分散は \(n \sigma^2\)であることはわかる。そのため、標準化した \[\begin{equation} \frac{S_n - n\mu}{\sqrt{n}\sigma} = \frac{\frac{1}{n}S_n - \mu}{\frac{\sigma}{\sqrt{n}}} \end{equation}\] は、平均0、分散1となることはわかる。しかしながら、この標準化した \(S_n\) が従う確率分布は未知である(はずである)。

しかしながら、中心極限定理では \[\begin{equation} \lim_{n \to \infty} P(a \le \frac{\frac{1}{n}S_n - \mu}{\frac{\sigma}{\sqrt{n}}} \le b) = \frac{1}{\sqrt{2\pi}} \int_a^b \exp \left(- \frac{u^2}{2} \right)du \end{equation}\] を保証する。すなわち、データ数が十分多いとき、\(S_n = X_1 + X_2 + ... + X_n\) を標準化したものはガウス分布に近づく。
データ数が十分多いという表現は曖昧ではあるが、便宜的にこの講義では \(N \ge 30\) のとき、データ数が十分多いとする。

この中心極限定理は、統計学の中でも最も重要な定理といっても過言ではない。なぜならば、平均、分散がある程度わかれば、どのようなデータでもガウス分布とみなせる変形が存在するためである。つまり、母集団分布は未知のままでも、変数変換した先はガウス分布とわかる。そのため、ガウス分布の性質さえ抑えておけば、母集団分布の分布自体はわからずとも問題ない。

中心極限定理をまず数値的に確かめてみよう。

dat = runif(10000)
hist(dat)

上記のように、0から1の値をとる一様分布に従う確率変数を100個生成してみる(平均0.5、分散1/12)。ついで中心極限定理と同様の変換を行う。これを1000回繰り返したヒストグラムと、平均0、分散1のガウス分布の確率分布を同時に描いてみる。

S = matrix(0, 1, 1000)
for(i in 1:1000){
  dat = runif(100)
  S[i] = (sum(dat)/length(dat) - 0.5)/(sd(dat)/sqrt(length(dat)))
}

hist(S, prob=T, xlim = c(-5, 5), ylim = c(0, 0.4))
curve(dnorm, -5, 5, add = T) # 確率密度関数

続いては、ガンマ分布と呼ばれる分布から生成してみる。

dat = rgamma(10000, 1, 1)
hist(dat)

このガンマ分布の平均、分散ともに1のため、

S = matrix(0, 1, 1000)
for(i in 1:1000){
  dat = rgamma(100, 1, 1)
  S[i] = (sum(dat)/length(dat) - 1)/(sd(dat)/sqrt(length(dat)))
}

hist(S, prob=T, xlim = c(-5, 5), ylim = c(0, 0.4))
curve(dnorm, -5, 5, add = T) # 確率密度関数

となる。確かに、中心極限定理が成り立っていることがわかる。つまり、データが多数計測可能である場合、標準化したデータ平均はガウス分布に従う。

中心極限定理の証明

\(\frac{S_n - n\mu}{\sqrt{n}\sigma}\)の分布がどのようなものになるかを計算する。分布を調べる手段として有効な指標が、モーメント母関数である。なぜならば、分散、歪度、尖度、など大多数の統計量を同時に検証できるためである。 \[\begin{equation} E \left[ \exp \left( t \frac{S_n - n\mu}{\sqrt{n}\sigma} \right)\right] = E \left[\exp \left( t \frac{X_1-\mu}{\sqrt{n}\sigma} \right)\exp \left( t \frac{X_2-\mu}{\sqrt{n}\sigma} \right) ... \exp \left( t \frac{X_n-\mu}{\sqrt{n}\sigma} \right)\right] \end{equation}\] ここで、すべてのデータが独立あること、そしてすべてのデータが同じ性質を有していることを利用すると、 \[\begin{align} E \left[ \exp \left( t \frac{S_n - n\mu}{\sqrt{n}\sigma} \right)\right] &= E \left[\exp \left( t \frac{X_1-\mu}{\sqrt{n}\sigma} \right)\right]E\left[\exp \left( t \frac{X_2-\mu}{\sqrt{n}\sigma} \right)\right] ... E\left[\exp \left( t \frac{X_n-\mu}{\sqrt{n}\sigma} \right)\right] \nonumber\\ &= \left( E \left[\exp \left( t \frac{X-\mu}{\sqrt{n}\sigma} \right)\right] \right)^n \end{align}\] となる。そして、最終的に補助変数 \(t\) を \(0\) に極限とることから \(t\) は十分小さいとして2次でテーラー展開を止めると、 \[\begin{align} \left( E \left[\exp \left( t \frac{X-\mu}{\sqrt{n}\sigma} \right)\right] \right)^n &\simeq \left( E \left[1 + t \frac{X-\mu}{\sqrt{n}\sigma} + \frac{t^2}{2!}\frac{(X-\mu)^2}{n\sigma^2} \right] \right)^n \nonumber\\ &= \left( 1 + \frac{t^2}{2n} \right)^n \end{align}\] ここで、 ネイピア数の定義 \(\mathrm{e} = \lim_{n\to\infty} \left( 1 + \frac{1}{n} \right)^n\) を思い出そう。すると、\(\frac{t^2}{2n} = \frac{1}{\hat{n}}\)とすると、 \[\begin{align} \lim_{n\to\infty} E \left[ \exp \left( t \frac{S_n - n\mu}{\sqrt{n}\sigma} \right)\right] &= \lim_{n\to\infty}\left( 1 + \frac{t^2}{2n} \right)^n \nonumber\\ &= \left( \lim_{\hat{n}\to\infty}\left( 1 + \frac{1}{\hat{n}} \right)^{\hat{n}} \right)^{\frac{t^2}{2}} \nonumber\\ &= \exp \left(\frac{t^2}{2} \right) \end{align}\]

さらに、ガウス分布のモーメント母関数が \(\exp\left(\frac{\sigma^2}{2}t^2 \right)\) であったことから、\(\frac{S_n - n\mu}{\sqrt{n}\sigma}\)は、\(\sigma^2 = 1\) としてガウス分布のモーメント母関数に \(N \to \infty\) にて収束することがわかる。つまり、分散 \(1\) のガウス分布になる。加えて、\(\frac{S_n - n\mu}{\sqrt{n}\sigma}\) の期待値は 0 になることはわかる \(E[S_n] = n\mu\)。

以上より、中心極限定理を証明できた。すなわち、\(\frac{S_n - n\mu}{\sqrt{n}\sigma}\)は、\(N \to \infty\) において、平均0、分散1のガウス分布に収束する。\(\frac{S_n - n\mu}{\sqrt{n}\sigma}\)は標準化そのもののため、平均0、分散1であることはわかる。中心極限定理はさらに、\(N \to \infty\) にて、標準化したデータ平均がガウス分布になることを保証してくれるものである。

演習問題