尤度のアニメーション

· 626 words · 2 minute read
R · bayes · psychology 
library(tidyverse)
## ─ Attaching packages ────────────────────────────────────────────────── tidyverse 1.3.0 ─
## ✓ ggplot2 3.3.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.3     ✓ dplyr   1.0.1
## ✓ tidyr   1.1.1     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ─ Conflicts ─────────────────────────────────────────────────── tidyverse_conflicts() ─
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

３つのデータポイントでの尤度

muが170だとすると？

datset <- c(165,173,182)
ggplot(data=data.frame(X=c(100,200)),aes(x=X))+
  stat_function(fun=dnorm,args = list(mean=170,sd=10),color="red")+
  geom_segment(aes(x = datset[1], xend = datset[1], y = 0, yend = dnorm(datset[1], 170, 10)), lty = 2) +
  geom_segment(aes(x = datset[2], xend = datset[2], y = 0, yend = dnorm(datset[2], 170, 10)), lty = 2) +
  geom_segment(aes(x = datset[3], xend = datset[3], y = 0, yend = dnorm(datset[3], 170, 10)), lty = 2) +
  geom_point(data=data.frame(X=datset,Y=rep(0,3)),aes(x=X,y=Y),size=3,color="blue") + 
  annotate("text", x=162,   y=0.001, label="165cm")+
  annotate("text", x=170,   y=0.001, label="173cm")+
  annotate("text", x=180,   y=0.001, label="182cm")

muが190だとすると？

ggplot(data=data.frame(X=c(100,200)),aes(x=X))+
  stat_function(fun=dnorm,args = list(mean=190,sd=10),color="red")+
  geom_segment(aes(x = datset[1], xend = datset[1], y = 0, yend = dnorm(datset[1], 190, 10)), lty = 2) +
  geom_segment(aes(x = datset[2], xend = datset[2], y = 0, yend = dnorm(datset[2], 190, 10)), lty = 2) +
  geom_segment(aes(x = datset[3], xend = datset[3], y = 0, yend = dnorm(datset[3], 190, 10)), lty = 2) +
  geom_point(data=data.frame(X=datset,Y=rep(0,3)),aes(x=X,y=Y),size=3,color="blue")+
  annotate("text", x=162,   y=0.001, label="165cm")+
  annotate("text", x=170,   y=0.001, label="173cm")+
  annotate("text", x=180,   y=0.001, label="182cm")

muが190だとすると？

ggplot(data=data.frame(X=c(100,200)),aes(x=X))+
  stat_function(fun=dnorm,args = list(mean=160,sd=10),color="red")+
  geom_segment(aes(x = datset[1], xend = datset[1], y = 0, yend = dnorm(datset[1], 160, 10)), lty = 2) +
  geom_segment(aes(x = datset[2], xend = datset[2], y = 0, yend = dnorm(datset[2], 160, 10)), lty = 2) +
  geom_segment(aes(x = datset[3], xend = datset[3], y = 0, yend = dnorm(datset[3], 160, 10)), lty = 2) + 
  geom_point(data=data.frame(X=datset,Y=rep(0,3)),aes(x=X,y=Y),size=3,color="blue")+
  annotate("text", x=162,   y=0.001, label="165cm")+
  annotate("text", x=170,   y=0.001, label="173cm")+
  annotate("text", x=180,   y=0.001, label="182cm")

平均値を動かすのは平行移動

正規分布が平行移動するのアニメーションを作るコードは次の通り。

animationパッケージとimage magickを別途ご用意ください。

library(animation)

ani.options(interval=0.1)
ani.options(loop=0)
ani.options(ani.width=800,ani.height=400)

likelihood_movie <- function(datset,fromX,toX){
  for(i in seq(fromX,toX,1)){
    g <- ggplot(data=data.frame(X=c(fromX,toX)),aes(x=X))+
      stat_function(fun=dnorm,args = list(mean=i,sd=10),color="red")+
      geom_segment(aes(x = datset[1], xend = datset[1], y = 0, yend = dnorm(datset[1], i, 10)), lty = 2) +
      geom_segment(aes(x = datset[2], xend = datset[2], y = 0, yend = dnorm(datset[2], i, 10)), lty = 2) +
      geom_segment(aes(x = datset[3], xend = datset[3], y = 0, yend = dnorm(datset[3], i, 10)), lty = 2) + 
      geom_point(data=data.frame(X=datset,Y=rep(0,3)),aes(x=X,y=Y),size=3,color="blue")+
      annotate("text", x=162,   y=0.001, label="165cm")+
      annotate("text", x=170,   y=0.001, label="173cm")+
      annotate("text", x=180,   y=0.001, label="182cm")+
      annotate("text", x=110,   y=0.038, label=paste0(" f(165 |",i,"10)=",dnorm(165,i,10)))+
      annotate("text", x=110,   y=0.036, label=paste0(" f(173 |",i,"10)=",dnorm(173,i,10)))+
      annotate("text", x=110,   y=0.034, label=paste0(" f(182 |",i,"10)=",dnorm(182,i,10)))
    
    print(g)
  }
}

saveGIF(likelihood_movie(datset=c(165,173,182),fromX=100,toX=200), 
        movie.name="likelihood.gif",clean=T)
## Output at: likelihood.gif
## [1] TRUE

出来上がった図がこちら。

likelihoodAnimation

Share