# (C) Petteri Pajunen 2006 # # # gibbsmixture.R # # simulate a mixture model with two Normal mixtures # # library(mvtnorm) # true parameter values lambda1 <- 0.25 lambda2 <- 0.75 mu1 <- matrix(c(-2,-2),2,1) mu2 <- matrix(c(2,2),2,1) cat("Variance scale? (0.1)") cvar <- scan(n=1,quiet=TRUE) Sigma <- cvar*matrix(c(1,0,0,1),2,2) cat("Initial mu1? (0 1)") m1 <- scan(n=2,quiet=TRUE) cat("Initial mu2? (1 0)") m2 <- scan(n=2,quiet=TRUE) samples <- 40 x <- matrix(0,samples,2) # generate values in proportion to lambda's x1max <- round(samples*lambda1) x <- rmvnorm(samples,mean=mu2,sigma=Sigma) x[1:x1max,] <- rmvnorm(10,mean=mu1,sigma=Sigma) X11(width=5,height=5) plot(x[1:x1max,1],x[1:x1max,2],xlim=c(-3,3),ylim=c(-3,3),type="p",col="red") points(x[(x1max+1):samples,1],x[(x1max+1):samples,2],col="blue") points(mu1[1],mu1[2],lwd=3,col="red") points(mu2[1],mu2[2],lwd=3,col="blue") iter <- 20 # initial values for Gibbs #m1 <- c(0.1,-0.1) #m2 <- c(-0.1,0.1) l1 <- 0.5 l2 <- 0.5 lim <- matrix(runif(2*samples,min=-1.0,max=1.0),samples,2) lim <- (lim>0) lim[,2] <- 1-lim[,1] mu1 <- m1 mu2 <- m2 y <- x for (i in 1:iter) { # simulate Lim for (j in 1:samples) { p1 <- l1*dmvnorm(y[j,],mean=m1,sigma=Sigma) psum <- l1*dmvnorm(y[j,],mean=m1,sigma=Sigma)+l2*dmvnorm(y[j,],mean=m2,sigma=Sigma) p1 <- p1/psum # this is the probability that Li1=1 u <- runif(1) lim[j,] <- c(0,0) if (u