# Estimating a simple mixture of normals #install.packages("geoR") # We need the rinvchisq() function for simulating random draws from the scaled inv chi-square. library(geoR) ########## BEGIN USER INPUT ################# # Data options #rawData <- read.table('/home/mv/Dropbox/Teaching/BayesLearn2012/Code/CanadianWages.dat', header = TRUE) #x <- as.matrix(rawData['logWage']) rawData <- geyser x <- as.matrix(rawData['duration']) # Model options nComp <- 3 # Number of mixture components # Prior options alpha <- 10*rep(1,nComp) # Dirichlet(alpha) muPrior <- rep(0,nComp) tau2Prior <- rep(100,nComp) sigma2_0 <- rep(var(x),nComp) nu0 <- rep(10,nComp) # MCMC options nIter <- 1000 # Number of Gibbs sampling draws # Plotting options plotFit <- TRUE lineColors <- c("blue", "green", "magenta", 'yellow') ################ END USER INPUT ############### ####### Defining a function that simulates from a Dirichlet distribution rDirichlet <- function(param){ nCat <- length(param) thetaDraws <- matrix(NA,nCat,1) for (j in 1:nCat){ thetaDraws[j] <- rgamma(1,param[j],1) } thetaDraws = thetaDraws/sum(thetaDraws) # Diving every column of ThetaDraws by the sum of the elements in that column. return(thetaDraws) } # Simple function that converts between two different representations of the mixture allocation S2alloc <- function(S){ n <- dim(S)[1] alloc <- rep(0,n) for (i in 1:n){ alloc[i] <- which(S[i,] == 1) } return(alloc) } # Initial value for the MCMC nObs <- length(x) S <- t(rmultinom(nObs, size = 1 , prob = rep(1/nComp,nComp))) # nObs-by-nComp matrix with component allocations. theta <- quantile(x, probs = seq(0,1,length = nComp)) sigma2 <- rep(var(x),nComp) probObsInComp <- rep(NA, nComp) # Setting up the plot xGrid <- seq(min(x)-1*apply(x,2,sd),max(x)+1*apply(x,2,sd),length = 100) xGridMin <- min(xGrid) xGridMax <- max(xGrid) mixDensMean <- rep(0,length(xGrid)) effIterCount <- 0 ylim <- c(0,2*max(hist(x)\$intensities)) for (k in 1:nIter){ message(paste('Iteration number:',k)) alloc <- S2alloc(S) # Just a function that converts between different representations of the group allocations nAlloc <- colSums(S) print(nAlloc) # Update components probabilities w <- rDirichlet(alpha + nAlloc) # Update theta's for (j in 1:nComp){ precPrior <- 1/tau2Prior[j] precData <- nAlloc[j]/sigma2[j] precPost <- precPrior + precData wPrior <- precPrior/precPost muPost <- wPrior*muPrior + (1-wPrior)*mean(x[alloc == j]) tau2Post <- 1/precPost theta[j] <- rnorm(1, mean = muPost, sd = sqrt(tau2Post)) } # Update sigma2's for (j in 1:nComp){ sigma2[j] <- rinvchisq(1, df = nu0[j] + nAlloc[j], scale = (nu0[j]*sigma2_0[j] + sum((x[alloc == j] - theta[j])^2))/(nu0[j] + nAlloc[j])) } # Update allocation for (i in 1:nObs){ for (j in 1:nComp){ probObsInComp[j] <- w[j]*dnorm(x[i], mean = theta[j], sd = sqrt(sigma2[j])) } S[i,] <- t(rmultinom(1, size = 1 , prob = probObsInComp/sum(probObsInComp))) } # Printing the fitted density against data histogram if (plotFit && (k%%10 ==0)){ effIterCount <- effIterCount + 1 hist(x, breaks = 20, freq = FALSE, xlim = c(xGridMin,xGridMax), main = paste("Iteration number",k), ylim = ylim) mixDens <- rep(0,length(xGrid)) components <- c() for (j in 1:nComp){ compDens <- dnorm(xGrid,theta[j],sd = sqrt(sigma2[j])) mixDens <- mixDens + w[j]*compDens lines(xGrid, compDens, type = "l", lwd = 2, col = lineColors[j]) components[j] <- paste("Component ",j) } mixDensMean <- ((effIterCount-1)*mixDensMean + mixDens)/effIterCount lines(xGrid, mixDens, type = "l", lwd = 3, col = 'red') legend("topleft", box.lty = 1, legend = c("Data histogram",components, 'Mixture'), col = c("black",lineColors[1:nComp], 'red'), lwd = 2) Sys.sleep(0) } } hist(x, breaks = 20, freq = FALSE, xlim = c(xGridMin,xGridMax), main = "Final fitted density") lines(xGrid, mixDensMean, type = "l", lwd = 2, col = "red") lines(xGrid, dnorm(xGrid, mean = mean(x), sd = apply(x,2,sd)), type = "l", lwd = 2, col = "blue") legend("topright", box.lty = 1, legend = c("Data histogram","Mixture density","Normal density"), col=c("black","red","blue"), lwd = 2)