廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

目で見る尤度関数(『ベイズ統計の理論と方法』より)

今日の川柳

ベイズ統計の理論と方法』1.4節の例を R でやってみます。

ベイズ統計の理論と方法

ベイズ統計の理論と方法

尤度関数が正規分布で近似できるとき、いろいろ良い性質がなりたちます。

関数 \phi(x) を標準正規分布の密度関数とし、確率モデル

 p(x)=(1-a) \phi(x) + a \phi(x-b)
(ただし  0 < a< 1) から100個の乱数を生成して、尤度関数の形を見ていきます。

真のパラメータが、(a_0,b_0)=(0.5,3) のとき、ヒストグラムを書いてみると、密度関数は混合分布で表せそうに見えます。

f:id:abrahamcow:20180330235444p:plain

最尤推定量は (0.52, 2.90) であり、真のパラメータの近くにあります。

尤度関数の等高線に二変量正規分布から生成した乱数を重ねてみると、尤度関数は正規分布で近似できそうに見えます。

f:id:abrahamcow:20180330235753p:plain

真のパラメータが、(a_0,b_0)=(0.5,1) のとき、ヒストグラムを書いてみると、ぱっと見では混合分布だとはわからないかもしれません。

f:id:abrahamcow:20180331000808p:plain

最尤推定量は (0.45 1.17) であり、まあまあ真のパラメータの近くにあります。

尤度関数の等高線はちょっといびつな形をしており、二変量正規分布で近似するのは少々強引な気がします。

f:id:abrahamcow:20180331001027p:plain

真のパラメータが、(a_0,b_0)=(0.5,0.5) のとき、ヒストグラムを書いてみると、やはりぱっと見では混合分布だとはわからないかもしれません。

f:id:abrahamcow:20180331001209p:plain

最尤推定量は (0.31 0.67) でした。真のパラメータからは少しずれています。

尤度関数の等高線は局所的にするどいピークを持ち、正規分布で近似するのはあきらかに無理があります。

f:id:abrahamcow:20180331001413p:plain

以下に R のコードを貼ります。

library(mvtnorm)
library(tidyverse)
rmixnorm <- function(n,a,b){
  Y <-sample.int(2, n, replace=TRUE, prob=c(1-a,a))
  mu <- c(0, b)
  rnorm(n, mu)
}
dmixnorm <- function(par,x,log=TRUE){
  a0 <- par[1]
  b0 <- par[2]
  ll <- sum(log((1-a0)*dnorm(x)+a0*dnorm(x,b0)))
  if(log){
    ll
  }else{
    exp(ll)
  }
}
set.seed(1)
x1 <- rmixnorm(100,0.5,3.0)
x2 <- rmixnorm(100,0.5,1.0)
x3 <- rmixnorm(100,0.5,0.5)
hist(x1)
hist(x2)
hist(x3)
opt1 <- optim(c(0.5,1), dmixnorm,control = list(fnscale=-1), x=x1, method = "L-BFGS-B",
             lower=c(0,-5),upper=c(1,5),hessian = TRUE)
opt2 <- optim(c(0.5,1), dmixnorm,control = list(fnscale=-1), x=x2, method = "L-BFGS-B",
             lower=c(0,-5),upper=c(1,5),hessian = TRUE)
opt3 <- optim(c(0.5,1), dmixnorm,control = list(fnscale=-1), x=x3)
print(round(opt1$par,2))
print(round(opt2$par,2))
print(round(opt3$par,2))
a <- seq(0,1,length.out = 200)
b <- seq(-5,5,length.out = 200)
parms <-expand.grid(a,b)
L1 <- apply(parms,1,dmixnorm,x=x1,log=FALSE)
datL1 <- as_data_frame(parms) %>% 
  set_names(c("a","b")) %>% 
  mutate(L=L1)
sample_out1 <-as_data_frame(rmvnorm(10000,opt1$par,-solve(opt1$hessian))) %>% 
  set_names(c("a","b"))
ggplot(datL1,aes(x=a,y=b))+
  geom_point(data=sample_out1,alpha=0.2)+
  geom_contour(aes(z=L),size=1)+
  theme_bw()

L2 <- apply(parms,1,dmixnorm,x=x2,log=FALSE)
datL2 <- as_data_frame(parms) %>% 
  set_names(c("a","b")) %>% 
  mutate(L=L2)
sample_out2 <-as_data_frame(rmvnorm(10000,opt2$par,-solve(opt2$hessian))) %>% 
  set_names(c("a","b"))
ggplot(datL2,aes(x=a,y=b))+
  geom_point(data=sample_out2,alpha=0.2)+
  geom_contour(aes(z=L),size=1)+
  theme_bw()

L3 <- apply(parms,1,dmixnorm,x=x3,log=FALSE)
datL3 <- as_data_frame(parms) %>% 
  set_names(c("a","b")) %>% 
  mutate(L=L3)
ggplot(datL3,aes(x=a,y=b))+
  geom_contour(aes(z=L),size=1)+
  theme_bw()




目で見る尤度関数