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()
3つのデータポイントでの尤度
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