先週のおさらい

(先々週のおさらい)連続な確率変数\(X\)の実現値が\((x, x+\Delta x) \ (\Delta x \ll 1)\)の範囲に収まる確率は\(p(X=x)\Delta x\)で定義され、\(p(X = x)\)を省略して\(p(x)\)と書く。\(p(X=x)\)の正式名称は確率密度関数であるが、確率分布と呼ばれることも多い。





これまでデータ\(\boldsymbol{x} = (x_1, ..., x_N)\)の平均は\(\bar{x} = \frac{1}{N}\sum_{i=1}^Nx_i\)として計算すること学んできた。一方、確率密度関数\(p(X = x)\)に従う確率変数の平均は\(E[X] = \bar[x] = \int_D xp(x)dx\)(連続)もしくは\(E[X] = \bar{x} = \sum_D xp(x)\)(離散)であることを学んだ。データ平均と確率密度関数の平均は一見異なるように見えるが、例えばコイントスを思い浮かべれば一致することがわかる。要は、\(x_i\)は\(p(x)\)に従い実現するため、両者は\(N\to\infty\)で完全に一致し、\(N\)が有限であれば誤差を伴い期待値として一致する(後々解説する)。





以下のガウス分布を学んだ。 \[\begin{equation} p(x | \mu, \sigma) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp \left( - \frac{1}{2\sigma^2}(x - \mu)^2 \right) \end{equation}\] 平均が\(\mu\), 分散が\(\sigma^2\)であることを計算した。さらに、歪度、尖度など3乗、4乗、、、の期待値もまた確率分布の代表値になることを学んだ。平均、分散、歪度、尖度をまとめて計算できる枠組みとして、モーメント母関数 \[\begin{equation} E[\exp(tX)] \end{equation}\] を学んだ。





統計のモチベーション

  統計学のモチベーションの一つが、計測データから母集団全体の特徴を推測することにある。 概念図もしくは対応表は下記の通り。













  重要な点は、
* 未知の計測不可能な母集団分布のデータは確率変数\(X\)であると想定する。
* 一部のデータを満遍なく(ランダムに)計測して、実現値\(x\)を得る。実現値は確率分布\(p(X = x)\)に従う。
* 計測したデータに基づき、未知の母集団全体の特徴を推測する。
ことである。特に、母集団分布の平均 \(\mu\) (未知)、母集団分布の分散 \(\sigma^2\) (未知) を推測することが多い。

  この表から、母集団の値は計測される前の確率変数\(X\)、データは計測された後の実現値\(x\)とみなすことができる。つまり、統計では、母集団の性質[平均 \(\mu\) (未知)、分散 \(\sigma^2\) (未知)]をもつ確率変数\(X\)を計測すると想定して、確率変数の実現値\(X = x\)を計測したデータとみなして様々な推論を行う。

  加えて、統計手法は適切に利用されなくてはいけない。母集団全体が計測できるならば計測する方が良い。満遍なく計測できていないデータはそれ相応の処理が必要である。

母集団分布

  上記のように、統計学では未知の母集団全体からランダムに計測したデータ(標本)を解析する。ランダムに計測する以上、計測する前のデータ \(X\) は確率変数とみなすことができ、その確率変数が従う確率密度関数 \(p(X = x)\) を母集団分布と呼ぶ(\(x\)は実現値もしくは計測されたデータ)。ランダムに計測する以上、\(i\)回目に計測されるであろうデータ \(X_i\) と \(j\) 回目に計測されるであろうデータ\(X_j\)は互いに独立であり、同じ \(p(X = x)\) からサンプリングされたものとみなすことができる(以下、\(p(X = x)\)を\(p(x)\)と略すこともある)。





  すなわち、確率変数の定義域 \(D\) において、\(\int_D x p(x) dx = \mu\)、\(\int_D (x - \mu)^2 p(x) dx = \sigma^2\)のとき (平均、分散だけ指定して、ほかは特に指定しない確率密度関数 \(p(x)\) を想定する)、
\[\begin{align} \int_D\int_D x_i x_j p(x_i, x_j) dx_i dx_j &= \int_D\int_D x_i x_j p(x_i)p(x_j) dx_i dx_j \ (\mathrm{独立})\nonumber\\ &= \int_D x_ip(x_i)dx_i \int_D x_jp(x_j) dx_j = \mu^2 \end{align}\] を満たす。これは重要な式なので、ぜひ抑えておいて欲しい。
以降、簡単のため、\(\int_D xp(x)dx = E[x]\)、\(\int_D x^2p(x)dx = E[x^2]\)、…として、\(E[\cdot]\) という記号を利用して期待値を表記することが多い。





変数変換と平均、分散

  \(a, b \in \mathbb{R}\) にて、確率変数 \(X\) を、\(aX + b\) と変換した際の平均、分散を計算してみよう。
\[\begin{equation} E[aX + b] = a \int_D x p(x)dx + b = a\mu + b \end{equation}\]





\[\begin{equation} E[(ax + b - [a\mu + b])^2] = a^2 E[(x - \mu)^2] = a^2 \sigma^2 \end{equation}\]





したがって、確率変数 \(X\) を、\(aX + b\) と変換した際の平均、分散は各々 \(a\mu + b\)、 \(a^2 \sigma^2\) となる。
分散に関しては、\(\mathrm{Var}(x) = \sigma^2\)や \[\begin{equation} \mathrm{Var}(ax + b) = a^2 \sigma^2 \end{equation}\] と表記することにする。

  これらを確認してみよう。

mean(rnorm(10000,0,1)) # μ = 1, σ = 1 のガウス分布から10000個確率変数を生成 → 平均
## [1] -0.003924157

\(a = 2\)、\(b = 1\)とすると、

a = 2
b = 1
mean(a*rnorm(10000,0,1) + b) # μ = 1, σ = 1 のガウス分布から10000個確率変数を生成 → 平均
## [1] 0.9963461

\(a = 2\)、\(b = 0\)とすると、

a = 2
b = 0
mean(a*rnorm(10000,0,1) + b) # μ = 1, σ = 1 のガウス分布から10000個確率変数を生成 → 平均
## [1] -0.04933121

ヒストグラムを書いてみると、

library(MASS)

h1 = rnorm(10000,0,1)

a = 2
b = 0
h2 = (a*rnorm(10000,0,1) + b) # μ = 1, σ = 1 のガウス分布から10000個確率変数を生成 → 平均

a = 2
b = 1
h3 = (a*rnorm(10000,0,1) + b) # μ = 1, σ = 1 のガウス分布から10000個確率変数を生成 → 平均

truehist(h1, xlim = c(-8, 8), ylim = c(0, 1), prob = T, col="#ff000080")
par(new=T)
truehist(h2, xlim = c(-8, 8), ylim = c(0, 1), prob = T, col="#00ff0080")
par(new=T)
truehist(h3, xlim = c(-8, 8), ylim = c(0, 1), prob = T, col="#0000ff80")

ここで、赤色は\(a = 1\)、\(b = 0\)の場合、緑色は\(a = 2\)、\(b = 0\)、青色は\(a = 2\)、\(b = 1\)である。感覚はつかめるであろう。

母集団分布

  \(K\)人の人が\(N\)個のデータを計測する状況を想定する。計測されるであろうデータが確率変数という想定である以上、Aさんが計測するデータの平均、Bさんが計測するデータの平均、Cさんの平均、Dさんの平均、…、これらは全て異なることが多い。すなわち、\(K\)人分のデータの平均を利用して「平均値の分布」もしくは「データ平均の分布」を描くことができ、このデータ平均の分布や、データ分散の分布のことを標本分布と呼ぶ。






例えば、平均0、分散1のガウス分布を母集団と仮定して(本来、平均、分散、確率分布の種類は未知)シミュレーションしてみる。1000人が10個のデータを計測して、データ平均、データ分散の分布を描いてみると以下の通り。

library(MASS)

N = 10 # 10サンプル
Ite = 1000 # 1000回繰り返す
  
mu = matrix(0, 1, Ite)
sigma = matrix(0, 1, Ite)
for(i in 1:Ite){
  dat = rnorm(N, 0, 1)
  mu[i] = mean(dat)
  sigma[i] = sd(dat)
}

par(mfrow=c(1,2)) 
truehist(mu, col="#ff000080")
truehist(sigma, col="#00ff0080")

やはり平均、分散の値は毎回異なり、標本分布を構成している。つまり、計測される前は、データ平均もデータ分散も標本分布に従う確率変数であることがわかる。計測のたびに、もしくは計測者によって異なる値を得られるデータ平均、データ分散からどのようにして高精度な推測値を得るか、というのが後半の講義のテーマである。

データ数と標本分布

  いくつデータを計測するかという問題は非常に根深い問題である。ある程度適切にデータ数を決定する方法はあるものの(power analysisなど)、この講義ではあまり深く踏み込まない。

ここでは、データ数とデータの平均、分散の関係性を検証してみる。以下は、1000人の人が1個のデータを計測した場合(赤)、5個のデータを計測した場合(緑)、100個のデータを計測した場合(黒)を想定している。そして各計測者がデータ平均を計算して得られた1000個のデータ平均のヒストグラムを表示している。もしくは、同一の計測者が1000回の計測を繰り返したとみなしてもよい。

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

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

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

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

  mu4[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="#00000020", add = T)

データ数が多いときほど、真の値0 (本当は未知の値) に近いことがわかる。

大数の法則

  上記のように、データ数が多いときほど、標本分布(もしくはデータ平均の分布)の平均は真の値に近づくことがわかる。






変数変換と平均値、分散の関係性から、データ平均の平均、分散を計算すると、 \[\begin{equation} E[\frac{1}{N}\sum_{i=1}^N x_i] = \frac{1}{N}\sum_{i=1}^N E[x_i] = \frac{1}{N}\sum_{i=1}^N \mu = \mu \end{equation}\] \[\begin{equation} \mathrm{Var}[\frac{1}{N}\sum_{i=1}^N x_i] = \frac{1}{N^2}\sum_{i=1}^N\mathrm{Var}[x_i] = \frac{1}{N}\sigma^2 \end{equation}\] すなわち、データ平均の期待値母集団の真の平均 \(\mu\) と一致する。加えて、データ平均の分散は \(\frac{1}{N}\sigma^2\) としてデータ数 \(N\) が増えると小さくなる。この傾向は確かに、上記のヒストグラムから確認できる。






  特に、\(N \to \infty\) にて、データ平均の分散が 0 になること、すなわち確率1にてデータ平均が母集団の真の平均に一致することを、大数の法則と呼ぶ。





不偏推定量

  上記の計算から、データ平均の平均は、 \[\begin{equation} E[\frac{1}{N}\sum_{i=1}^N x_i] = \mu \end{equation}\] となる。すなわち、計測したデータの平均は、期待値的に本来未知の母集団の真の平均と一致する。言い換えるならば、計測したデータから平均を計算してこれを未知の母集団の平均の推定値とすると、期待値としては本来未知の母集団の真の平均と一致する。データ平均は真の未知の平均の推定値としては“良い”ものであることがわかる。






  このように、期待値が真の母集団の特徴と一致する推定量のことを、不偏推定量と呼ぶ。データの平均は、未知の母集団の平均の不偏推定量である。ちなみに、推定量とは平均や分散などの用語そのものと同等、推定値は0.1や0.005など具体的な数値を意味する。

標準誤差

  \[\begin{equation} \mathrm{Var}[\frac{1}{N}\sum_{i=1}^N x_i] = \frac{1}{N^2}\sum_{i=1}^N\mathrm{Var}[x_i] = \frac{1}{N}\sigma^2 \end{equation}\] のルートをとったものである \[\begin{equation} \frac{\sigma}{\sqrt{N}} \end{equation}\] を、標準誤差と呼ぶ。標準誤差は、データ平均のばらつきを意味し、データ平均を真の母集団平均の推定値として採用する上で、データ平均と真の母集団平均の間にどれほどの違いが有り得るかを示す。データ平均と直接足し算引き算をして意味がある値は、この標準誤差である。論文においてデータ平均をプロットするときには、真の母集団平均の推定値を示すことが多く、(ほぼ必ず)標準誤差を足し引きしたものを同時にプロットする必要がある。以下は、100個のデータを計測→データ平均と標準誤差を計算、という流れを10回行いプロットした図。

N4 = 100
Ite = 10 # 1000回繰り返す
  
mu4 = matrix(0, 1, Ite)
sem4 = matrix(0, 1, Ite)

for(i in 1:Ite){
  dat = rnorm(N4, 0, 1)
  mu4[i] = mean(dat)
  sem4[i] = sd(dat)/sqrt(N4)
}

b = barplot(mu4, ylim = c(-0.35, 0.35)) 
arrows(b, c(mu4 - sem4), b, c(mu4 + sem4),length=0.05,angle=90, code = 3)

統計の概念図

  以上をまとめると、下記の表のとおりになる。











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





演習問題