演習問題

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





例えば、期末テストなどの場合(母集団のデータ全てを計測可能な場合)、クラス内のばらつきを知りたければ標準偏差で十分。センター試験の全国平均を1000人の得点から推測するとき(母集団の一部を観測した場合)、推定値はデータ平均\(\frac{1}{1000}\sum_{i = 1}^{1000}x_i\)が不偏推定量になり、推定値のばらつき(もしくは真の値からのおおよそのズレ)は標準誤差で表す。





標準誤差は \[\begin{equation} \mathrm{Var}\left(\frac{1}{N}\sum_{i = 1}^NX_i\right) \end{equation}\] であり、分散の変数変換に従えば\(\mathrm{Var}\left(\frac{1}{N}\sum_{i = 1}^NX_i\right) = \frac{\sigma^2}{N}\)。地道に計算すると、\(E[(X_i - \mu)(X_j - \mu)]\)より\((i \neq j)\)、 \[\begin{eqnarray} \mathrm{Var}\left(\frac{1}{N}\sum_{i = 1}^NX_i\right) &=& E \left[ \left( \sum_{i = 1}\frac{1}{N}X_i - \mu \right)^2 \right] = E \left[ \left( \frac{X_1}{N} + \frac{X_2}{N} + ... + \frac{X_N}{N} - \frac{\mu}{N} - \frac{\mu}{N} ... - \frac{\mu}{N} \right)^2 \right] \nonumber\\ &=& \frac{1}{N^2}E \left[ \sum_{i = 1}^N\sum_{j = 1}^N(X_i - \mu)(X_j - \mu) \right] = \frac{1}{N}E \left[\frac{1}{N}\sum_{i = 1}^N(X_i - \mu)^2 \right] = \frac{\sigma^2}{N} \end{eqnarray}\] データ平均\(\frac{1}{N}\sum_{i = 1}^NX_i\)と真の母集団平均\(\mu\)との２乗誤差を表していると言える。

先週のおさらい

  統計の想定している状況を学んだ。













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





計測する前のデータ(確率変数)は、母集団分布に従うと仮定する。未知の真の平均\(\mu\), 未知の真の分散\(\sigma^2\)を仮定して、確率変数は以下の性質に従うものとする。\(E[X_i] = \mu\), \(\mathrm{Var}[X_i] = \sigma^2\), \(E[X_iX_j] = 0\)。





確率変数 \(X\) を \(aX + b\) として変数変換すると、 \[\begin{equation} E[aX + b] = aE[X] + b = a\mu + b \end{equation}\] \[\begin{equation} \mathrm{Var}[aX + b] = a^2\mathrm{Var}[[X] = a^2\sigma^2 \end{equation}\] となることを学んだ。





これを利用すると、 \[\begin{equation} E[\frac{1}{N}\sum_{i = 1}^N] = \mu \end{equation}\] となり、データ平均の期待値は真の平均と一致すること(不偏推定量)、データ平均の分散は \[\begin{equation} \mathrm{Var}[\frac{1}{N}\sum_{i = 1}^N] = \frac{\sigma^2}{N} \end{equation}\] となりこれを標準誤差と呼ぶ。\(N \to \infty\)のとき、\(\frac{1}{N}\sum_{i = 1}^N \to \mu\)となり、これを大数の法則と呼ぶ。





標本分散 (データ分散)

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

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))

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

  上図は、赤色が\(N = 2\)、緑色が\(N = 5\)、青色が\(N = 10\)、黒色が\(N = 10000\)のときにデータ分散を計算し、この操作を1000回繰り返した際の分散のヒストグラム (このような分布を標本分布と呼ぶ) である。真の母集団分散は1 (シミュレーションなのでわかる) に対して、\(N = 10000\)のときはおおよそ1に一致している。 しかしながら、\(N = 2, 5, 10\)のときは、明らかに1とは異なることがわかる。






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





\[\begin{align} E\left[ \frac{1}{N}\sum_i \left( X_i - \frac{1}{N} \sum_j X_j \right)^2 \right] &= E\left[ \frac{1}{N}\sum_i \left( X_i - \mu + \mu - \frac{1}{N} \sum_j X_j \right)^2 \right] \nonumber\\ &= E\left[ \frac{1}{N}\sum_i \left( X_i - \mu \right)^2\right] + E\left[ \frac{1}{N}\sum_i \left(\mu - \frac{1}{N} \sum_j X_j\right)^2\right] + E\left[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[(X_i - \mu)(X_j - \mu)] = 0 \ (i \neq j)\)を利用すると、 \[\begin{equation} E\left[\frac{1}{N}\sum_i \left(\mu - \frac{1}{N} \sum_j X_j\right)^2\right] = E\left[\frac{1}{N}\sum_i \frac{1}{N^2}\sum_j\left(X_j - \mu\right)^2\right] = \frac{1}{N}E\left[\frac{1}{N}\sum_i \frac{1}{N}\sum_j\left(X_j - \mu\right)^2\right] = \frac{\sigma^2}{N} \end{equation}\] を得る。計算するべき値は \(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[X_iX_j] = 0 \ (i\neq j)\)、\(E[X_iX_i] = \sigma^2 + \mu^2\)に注意。











\[\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\) で割ることを実装していることが多いように見受けられる (そのため、分散の計算式を手打ちすると、ソフトウェアが出力した分散と異なるため困惑することがある)。

\(N-1\)で割るデータ分散の分布(標本分布)を計算してみる。

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))
t3 = c(mean(mu1)*N1/(N1-1), mean(mu2)*N2/(N2-1), mean(mu3)*N3/(N3-1), mean(mu4)*N4/(N4-1))

hist(mu1*N1/(N1-1), prob = T, ylim = c(0, 2), col="#ff000020")
hist(mu2*N2/(N2-1), prob = T, ylim = c(0, 2), col="#00ff0020", add = T)
hist(mu3*N3/(N3-1), prob = T, ylim = c(0, 1), col="#0000ff20", add = T)
#hist(mu4, prob = T, ylim = c(0, 2), col="#0000ff20", add = T)
hist(mu5*N5/(N5-1), prob = T, ylim = c(0, 2), col="#00000020", add = T)

赤色が\(N = 2\)、緑が\(N = 5\)、青が\(N = 10\)、黒が\(N = 10000\)の例を示している。標本分布ではわかりづらいため、標本分布の平均を下記に示す。

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))
t3 = c(mean(mu1)*N1/(N1-1), mean(mu2)*N2/(N2-1), mean(mu3)*N3/(N3-1), mean(mu4)*N4/(N4-1))

t1
## [1]   2   5  10 100
t2
## [1] 0.4777523 0.8147995 0.9110261 0.9945339
t3
## [1] 0.9555046 1.0184994 1.0122513 1.0045797

\(N = 2, 5, 10, 100\)において、上から\(N\)を用いた分散の標本分布の平均、その下は\(N-1\)を用いた分散の標本分布の平均。\(N-1\)を用いた場合には明らかに母集団分散\(1\)に近いことがわかる。





以上より、データ解析において、1) 真の母集団平均を推定したいならばデータ平均\(\bar{x}\)が適していること、2) 真の平均とデータ平均のずれは標準誤差\(\sigma^2/N\)で表すことができること、真の母集団分散\(\sigma^2\)の推定値は\(\frac{1}{N-1}\sum_{i = 1}^N(x_i - \bar{x})\)が適していることを明らかにした。
NOTE: これらはあくまで「不偏推定量」であり、推定方法により推定量は異なることに留意。実際、最尤推定という手法では、分散の推定値は\(\frac{1}{N}\sum_{i = 1}^N(x_i - \bar{x})\)。

標準化

  変数変換の知識を利用して、データ解析で最も重要な要素の一つである、前処理の方法を学ぶ。前処理は、データがいくつかあるときに、変動が大きいデータを重要視してしまうこと、大きな値をもつデータを重要視してしまうことを避けるために必要不可欠である。変動が大きいデータと変動が小さいデータを平等に扱うためには、標準誤差を1にすること、値の大きさを平等に扱うためには平均を0にすることが一つの手法である。標準誤差を1、平均を0にするための変数変換を求めればよい。





  その変換方法は、平均 \(\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\), … を公平に比較する際や同時に解析する際に利用する。





  (余談)しばしば受験にて利用される指標である偏差値は、この標準化を利用したもの。得点データの平均\(\mu\)、得点データの標準偏差\(\sigma\)を用いて、得点\(x\)のときの偏差値\(z\)は \[\begin{equation} z(x) = 50 + 10 \frac{x - \mu}{\sigma} \end{equation}\] となる。平均点と同じならば偏差値は50、そこから標準偏差だけ離れるごとに10ずつ変化する。





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

  確率変数\(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) # 確率密度関数

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

  続いては、二項分布から生成してみる。

dat = rbinom(10000, 10, 0.2)
hist(dat)

この二項分布の平均を0に合わせると、

S = matrix(0, 1, 1000)
for(i in 1:1000){
  dat = rbinom(10000, 10, 0.2)
  S[i] = (sum(dat)/length(dat) - 1)/(sd(dat)/sqrt(length(dat)))
}
mS = mean(S)
S = S - mS

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\) にて、標準化したデータ平均がガウス分布になることを保証してくれるものである。





演習問題