Skip to content

Instantly share code, notes, and snippets.

@abikoushi
Created December 8, 2023 16:37
Show Gist options
  • Save abikoushi/a0eaeecfd9e2d50c03cc96605c0e5bd4 to your computer and use it in GitHub Desktop.
Save abikoushi/a0eaeecfd9e2d50c03cc96605c0e5bd4 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")
@abikoushi
Copy link
Author

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment