posterior density of the mixture normal distribution

lpmixnorm.R

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")

log marginal likelihood can be calculate as follows:

logsumexpv <- function(x){
  mx <- max(x)
  mx + log(sum(exp(x-mx)))
}

logsumexp=function (logx1,logx2){
  logx1 + log1p(exp(logx2-logx1))
}
llmixnorm0 <- 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)))
}

N <- 100L
set.seed(1);y11 <- c(rnorm(N/2),rnorm(N/2,0.5))

theta=seq(0.5,3,by=0.1)
len <- length(theta)
lp <- numeric(len)
for(i in 1:len){
  Pmat <- cbind(rbeta(100000,theta[i],theta[i]),runif(100000))
  lp[i] <- logsumexpv(apply(Pmat, 1, llmixnorm0, y=y11))
}

png("logmarginal.png")
plot(theta,lp,type="l",pch=16,ylab="log-marginal likelihood")
dev.off()