Skip to content

Instantly share code, notes, and snippets.

@ysaito8015
Forked from abikoushi/lpmixnorm.R
Created December 26, 2023 02:17
Show Gist options
  • Save ysaito8015/eda8d7776db34f4fb18835952b9c8914 to your computer and use it in GitHub Desktop.
Save ysaito8015/eda8d7776db34f4fb18835952b9c8914 to your computer and use it in GitHub Desktop.
posterior density of the mixture normal distribution
library(ggplot2)
library(gganimate)
logsumexp=function (logx1,logx2){
logx1 + log1p(exp(logx2-logx1))
}
llmixnorm <- function(par, y){
a0 <- par[1]
b0 <- par[2]
theta <- par[3]
sum(logsumexp(log(1-a0)+dnorm(y,log = TRUE), log(a0)+dnorm(y,b0,log = TRUE))) + dbeta(a0, theta, theta, log = TRUE)
}
N <- 100L
set.seed(1)
y11 <- c(rnorm(N/2),rnorm(N/2,0.5))
ggplot()+
geom_histogram(data=NULL, aes(x=y11), fill="grey70", bins=25)+
theme_bw(14)+labs(x="y", y="count")
parms <- expand.grid(a=seq(0.3,0.99,length.out = 200),
b=seq(0,1,length.out = 200),
theta=seq(0.5,3,by=0.1))
l11 <- apply(parms, 1, llmixnorm, y=y11)
dfL11 <- data.frame(parms, value = exp(l11))
ggplot(dfL11,aes(x=b,y=a))+
geom_point(aes(colour=value))+
scale_colour_continuous(type = "viridis")+
labs(title = 'theta: {round(frame_time,2)}') +
transition_time(theta)+
theme_bw(14)+theme(legend.position = "none")
anim_save("post_mixnorm.gif")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment