正規分布のパラメタを推定する

  • ベルヌーイ乱数から平均と分散を推定することと、ベルヌーイ乱数からベルヌーイ分布のパラメタを(最尤)推定することが違っていたので
  • 正規乱数の場合にはどうなるかを見てみる
  • 正規分布は2変数分布なので、標本から母分布の平均と分散とを推定してやれば、それに合致した正規分布というのは、存在する。この点はベルヌーイ分布の推定の場合と違っている
# 正規乱数発生
m <- pi
s2 <- exp(1)
n <- 20
X <- rnorm(n,m,sqrt(s2))
# 平均と分散とでグリッド形成
ms <- seq(from=2,to=5,length=100)
s2s <- seq(from=0.5,to=6,length=100)
# 対数尤度を計算して行列状にする
LL <- matrix(0,length(ms),length(s2s))
for(i in 1:length(ms)){
	for(j in 1:length(s2s)){
		LL[i,j] <- sum(log(dnorm(X,ms[i],sqrt(s2s[j]))))
	}
}
image(ms,s2s,LL)
# 母分布の平均と分散は赤い点
points(m,s2,pch=20,col=2)
# 標本平均と標本からの不偏分散は緑の点
points(mean(X),var(X),pch=20,col=3)
# 標本平均と標本分散は青の点
points(mean(X),var(X)*(n-1)/n,pch=20,col=4)
contour(ms,s2s,LL,add=TRUE,nlevels=100)
sum(log(dnorm(X,mean(X),sqrt(var(X)))))
sum(log(dnorm(X,mean(X),sqrt(var(X)*(n-1)/n))))

  • 別の試行だが、拡大すると:標本分散が最尤推定になっているように見える

  • 対数尤度関数を解いて確かめてみる
  • LL(\mu,\sigma|X=(x_i)) = \log{\frac{1}{\sqrt{2\pi \sigma^2}}} -\frac{\sum_i (x_i-\mu)^2}{2\sigma^2}
  • \mu,\sigma偏微分すれば
  • \frac{\partial LL}{\partial \mu}= \frac{1}{2\sigma^2}\sum 2(x_i-\mu)=0から、\hat{\mu}=\frac{\sum_i x_i}{N}
  • \frac{\partial LL}{\partial \sigma} = -\frac{N}{\sigma} + \frac{\sum_i (x_i-\mu)^2}{\sigma^3} = \frac{1}{\sigma^3}(\sum_i(x_i-\mu)^2 -N \sigma^2) = 0から、\hat{\sigma^2} = \frac{1}{N}\sum_i (x_i-\mu)^2で、\hat{\mu}を採用すれば、\hat{\sigma^2}= \frac{1}{N} \sum_i (x_i-\hat{\mu})^2
  • 結局、標本平均を母分布の平均とし、標本分散を母分布の分散として推定するというのは、母分布が正規分布だと仮定した上での正規分布の平均パラメタと分散パラメタの最尤推定をしていることだとわかる
  • 逆に、不偏分散を分散の推定値とするということは、正規分布を仮定せずに推定している何か(分布についての事前仮定をしないときの最尤推定??)、ということになる
  • 対数尤度は標本分散を使う方が大きいが、座標「平均、分散」は標本からの不偏分散を用いる方が、この例では、母分布のそれに近い
> sum(log(dnorm(X,mean(X),sqrt(var(X)))))
[1] -32.21791
> sum(log(dnorm(X,mean(X),sqrt(var(X)*(n-1)/n))))
[1] -32.20498
  • 試行回数を増やしてみよう
    • 平均は変わらず、標本分散と不偏分散との関係は、常に標本分散の方が小だから、標本分散を赤、不偏分散を緑、母分布を黒にすると、こんな感じ

m <- pi
s2 <- exp(1)
n <- 20
n.iter <- 100

ms <- vars <- vars.ub <- rep(0,n.iter)
for(i in 1:n.iter){
	X <- rnorm(n,m,sqrt(s2))
	ms[i] <- mean(X)
	vars.ub[i] <- var(X)
	vars[i] <- vars.ub[i]*(n-1)/n
}

X <- c(m,ms,ms)
Y <- c(s2,vars,vars.ub)

plot(X,Y,col=c(1,rep(2,n.iter),rep(3,n.iter)),pch=20,cex=c(3,rep(1,n.iter*2)))
segments(ms,vars,ms,vars.ub)