となる確率分布を考える。この確率分布において
*確率分布の定義を満たしていること、
*\(x\)と\(y\)の平均、
*\(x\)と\(y\)の標準偏差、
*\(x\)と\(y\)の共分散、
*\(x\)と\(y\)の相関係数
を計算せよ。
*明らかに\(p(x, y) \ge 0\)である。加えて、
\(\int_0^1\int_0^1 p(x, y)dxdy = \int_0^1 \left[ \frac{x^2}{2} \right]_0^1+ y[x]_0^1) dy = \int_0^1 (\frac{1}{2} + y) dy = 1\)
したがって、確率分布の定義を満たす。
*\(x\)と\(y\)は入れ替え対称なので、平均も分散も標準偏差も同じ。まずは\(x\)の平均を計算する。
\(\int_0^1\int_0^1 xp(x, y)dxdy = \int_0^1 (\left[ \frac{x^3}{3} \right]_0^1+ y\left[\frac{x^2}{2} \right]_0^1)dy = \int_0^1 (\frac{1}{3} + \frac{y}{2}) dy = \frac{7}{12}\)
*ついで\(x^2\)の期待値を計算する。
\(\int_0^1\int_0^1 x^2p(x, y)dxdy = \int_0^1 (\left[ \frac{x^4}{4} \right]_0^1+ y\left[\frac{x^3}{3} \right]_0^1)dy = \int_0^1 (\frac{1}{4} + \frac{y}{3}) dy = \frac{5}{12}\)
したがって、分散は
\(\mathrm{Var}(x) = \frac{5}{12} - (\frac{7}{12})^2 = \frac{11}{144}\)
標準偏差は \(\mathrm{Std}(x) = \frac{\sqrt{11}}{12}\)
*共分散を計算するために、\(xy\)の期待値を計算する。 \(\int_0^1\int_0^1 xy\ p(x, y)dxdy = \int_0^1 (\left[ \frac{x^3}{3} \right]_0^1y+ y^2\left[\frac{x^2}{2} \right]_0^1)dy = \int_0^1 (\frac{y}{3} + \frac{y^2}{2}) dy = \frac{1}{3}\)
どちらも平均は\(\frac{7}{12}\)であり、共分散は \(\frac{1}{3} - \frac{49}{144} = -\frac{1}{144}\)
*相関係数は、 \(\frac{-\frac{1}{144}}{\frac{\sqrt{11}}{12}\frac{\sqrt{11}}{12}} = -\frac{1}{11}\)
連続な確率変数\(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}\] \(f(x)\)の\(p(x)\)に関する期待値を\(E[f(x)] = \int_D f(x)p(X)dx\)とすると(\(D\)は\(x\)の定義域)、ガウス分布の平均は \[\begin{equation}
E[x] = \frac{1}{\sqrt{2\pi\sigma^2}} \int_{-\infty}^{\infty} x \exp \left( - \frac{1}{2\sigma^2}(x - \mu)^2 \right)dx
\end{equation}\] として計算できる。この積分を計算すると、ガウス分布の平均は\(\mu\)であることがわかる。
これは分散の定義通り計算した方が見通しが良い。 \[\begin{equation}
\mathrm{Var}(x) = E[(x-\mu)^2] = \frac{1}{\sqrt{2\pi\sigma^2}} \int_{-\infty}^{\infty} (x-\mu)^2 \ \exp \left( - \frac{1}{2\sigma^2}(x-\mu)^2\right) dx
\end{equation}\] \(\hat{x} = x - \mu\) として、改めて \(\hat{x}\) を \(x\) と表記すると、 \[\begin{equation}
E[x^2] = \frac{1}{\sqrt{2\pi\sigma^2}} \int_{-\infty}^{\infty} x^2 \ \exp \left( - \frac{1}{2\sigma^2}x^2\right) dx
\end{equation}\]
\(x^2 \ \exp \left( - \frac{1}{2\sigma^2}x^2\right) = -\sigma^2 x \ \frac{d}{dx} \exp \left( - \frac{1}{2\sigma^2}x^2\right)\) であることに着目して、部分積分を利用すると (\((fg)' = f'g + fg'\) より、 \(\int fg' = [fg] - \int f'g\)) 、
\[\begin{align}
E[x^2] &= -\sigma^2 \left( [x \ \exp\left( - \frac{1}{2\sigma^2}x^2\right)]_{-\infty}^{\infty} - \frac{1}{\sqrt{2\pi\sigma^2}} \int_{-\infty}^{\infty} \exp \left( - \frac{1}{2\sigma^2}x^2\right) dx\right) \nonumber\\
&= \sigma^2
\end{align}\]
したがって、ガウス分布の分散は \(\sigma^2\) である。これも確かめてみよう。
var(rnorm(10000,1,1)) # μ = 1, σ = 1 のガウス分布から10000個確率変数を生成 → 分散
## [1] 1.011489
var(rnorm(10000,2,1)) # μ = 2, σ = 1 のガウス分布から10000個確率変数を生成 → 分散
## [1] 0.9939311
var(rnorm(10000,1,2)) # μ = 1, σ = 2 のガウス分布から10000個確率変数を生成 → 分散 (σ×σ)
## [1] 3.953751
以上より、たしかにガウス分布の分散は \(\sigma^2\) である。
これまでに、ベルヌーイ分布、二項分布、一様分布、ガウス分布など、様々な確率分布を学んできた。 多くの場合、平均、分散の2つを定量化することが多い。しかしながら、確率分布を表現する特徴量はこの2つだけではない。
その他の確率分布を特徴づける値として、歪度や尖度が挙げられる。歪度とは、平均 \(\mu\) において \(\frac{(x-\mu)^3}{\sigma^3}\) の期待値、尖度とは \(\frac{(x-\mu)4}{\sigma^4}\) を意味する。
歪度を計算してみると、下記の通り。
dam = rnorm(100000,1,1) # μ = 1, σ = 1 のガウス分布から10000個確率変数を生成
mean((dam - mean(dam))^3)
## [1] -0.0051974
dam = rnorm(100000,1,2) # μ = 1, σ = 2 のガウス分布から10000個確率変数を生成
mean((dam - mean(dam))^3)
## [1] 0.05548251
dam = rnorm(100000,2,1) # μ = 1, σ = 1 のガウス分布から10000個確率変数を生成
mean((dam - mean(dam))^3)
## [1] 0.001293371
上記のように、ガウス分布の歪度は0になる (解析的に0)。
尖度を計算してみると、下記の通り。
dam = rnorm(100000,1,1) # μ = 1, σ = 1 のガウス分布から10000個確率変数を生成
mean((dam - mean(dam))^4)
## [1] 2.984805
dam = rnorm(100000,1,2) # μ = 1, σ = 2 のガウス分布から10000個確率変数を生成
mean((dam - mean(dam))^4/((2^2)^2))
## [1] 2.987444
dam = rnorm(100000,2,1) # μ = 1, σ = 1 のガウス分布から10000個確率変数を生成
mean((dam - mean(dam))^4)
## [1] 3.026457
dam = rnorm(100000,2,3) # μ = 1, σ = 1 のガウス分布から10000個確率変数を生成
mean((dam - mean(dam))^4/((3^2)^2))
## [1] 2.957002
上記のように、ガウス分布の尖度は、\(3\)になる。
もちろんついで5乗, 6乗, … と計算していくことも可能ではあるが、どう考えても効率が悪い。 そこで確率分布を特徴づける方法として、モーメント母関数が提案された。
モーメント母関数とは、補助変数\(t\)を用いて
\[\begin{equation}
E[e^{t(X-\mu)}]
\end{equation}\] として定義される(正確には、これは中心化モーメント母関数とも呼ばれ、\(E[e^{t(X)}]\)をモーメント母関数と呼ぶこともある)。
このモーメント母関数をテーラー展開すると、
\(E[e^{t(X-\mu)}] = E \left[1 + t(X-\mu) + \frac{1}{2!}(t(X-\mu))^2 + \frac{1}{3!}(t(X-\mu))^3 + ... \right]\)
として無限級数で書き表すことができる。
そして、例えば \(t\) で微分して \(\lim_{t\to0}\) とすると、 \(\lim_{t\to0}\frac{d}{dt}E[e^{t(X-\mu)}] = E[(X-\mu)]\)
となる。最初に \(\mu = 0\) としてこの値を計算すると、平均 \(E[X]\) を得る。
ついで\(t^2\)で微分して \(\lim_{t\to0}\) とすると、
\(\lim_{t\to0}\frac{d^2}{dt^2}E[e^{t(X-\mu)}] = E[(X-\mu)^2]\)
つまり、分散となる。このようにして、 \(\lim_{t\to0}\frac{d^k}{dt^k}E[e^{tX}] = E[(X - \mu)^k]\)
を得る。つまりモーメント母関数さえ計算してしまえば、分散、歪度、尖度を始めとして、非常に多数の特徴量を計算可能になる。加えて、モーメント母関数が一致する2つの確率分布は、必ず等しい確率分布になるという性質ももつ。
ガウス分布のモーメント母関数を計算する。定義より、
\[\begin{equation}
E[e^{tX}] = \frac{1}{\sqrt{2\pi\sigma^2}} \int_{-\infty}^{\infty} \exp(t(x-\mu))\exp( - \frac{1}{2\sigma^2} (x-\mu)^2)dx
\end{equation}\] となる。指数関数の中身を整理すると、
\(tx - t\mu - \frac{1}{2\sigma^2}x^2 + \frac{\mu}{\sigma^2}x - \frac{1}{2\sigma^2}\mu^2 = - \frac{1}{2\sigma^2}x^2 + \left(t + \frac{\mu}{\sigma^2}\right)x - t\mu - \frac{1}{2\sigma^2}\mu^2\)
となる。ここで平方完成すると、
\(- \frac{1}{2\sigma^2}(x - \sigma^2 \left(t + \frac{\mu}{\sigma^2}\right) )^2 + \frac{\sigma^2}{2}\left( t^2 + 2\frac{\mu}{\sigma^2}t + \frac{\mu^2}{\sigma^4}\right) - t\mu - \frac{1}{2\sigma^2}\mu^2 = - \frac{1}{2\sigma^2}(x - \sigma^2 \left(t + \frac{\mu}{\sigma^2}\right) )^2 + \frac{\sigma^2}{2}t^2\)
を得る。つまり、
\[\begin{align}
E[e^{tX}] &= \frac{1}{\sqrt{2\pi\sigma^2}} \int_{-\infty}^{\infty} \exp\left(\frac{1}{2\sigma^2}(x - \sigma^2 \left(t + \frac{\mu}{\sigma^2}\right) )^2\right)\exp\left(\frac{\sigma^2}{2}t^2 \right)\nonumber\\
&= \exp\left(\frac{\sigma^2}{2}t^2 \right)dx
\end{align}\]
ここから分散、歪度、尖度…などを求めることができる。
統計学のモチベーションの一つが、計測データから母集団全体の特徴を推測することにある。 概念図は下記の通り。
重要な点は、
* 未知の計測不可能な母集団分布から、
* 一部のデータを満遍なく(ランダムに)計測し、
* 計測したデータに基づき、未知の母集団全体の特徴を推測する
ことである。特に、母集団分布の平均 \(\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.0174935
\(a = 2\)、\(b = 1\)とすると、
a = 2
b = 1
mean(a*rnorm(10000,0,1) + b) # μ = 1, σ = 1 のガウス分布から10000個確率変数を生成 → 平均
## [1] 1.014048
\(a = 2\)、\(b = 0\)とすると、
a = 2
b = 0
mean(a*rnorm(10000,0,1) + b) # μ = 1, σ = 1 のガウス分布から10000個確率変数を生成 → 平均
## [1] 0.03964506
ヒストグラムを書いてみると、
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\)も未知)。
データ \(x_1, x_2, ..., x_N\) を計測したとき (平均 \(\mu\) )、このデータの標準偏差 \(\sigma = \sqrt{\frac{1}{N} \sum_{i=1}^N(x_i - \mu)^2}\) と標準偏差 \(\frac{\sigma}{N}\) を両方計算することができる。標準偏差と標準誤差の違いは何かを述べよ。
ガウス分布のモーメント母関数を求めよ (次の講義で利用するので、前回の講義の計算をもう一度自分で確かめよう)