# Mixture of Gaussians on Log of Data

I am practicing Mixture of Gaussians and found the below dataset snoq, which is the precipitation amounts recorded at a US region, with NA and no precipitation days removed.

snoqualmie <- read.csv("http://www.stat.cmu.edu/~cshalizi/402/lectures/16-glm-practicals/snoqualmie.csv",header=FALSE)
snoqualmie.vector <- na.omit(unlist(snoqualmie)) # remove NA's and flatten
snoq <- snoqualmie.vector[snoqualmie.vector > 0] # days where precipitation was greater than 0


In the exercise code, the instructor fits a mixture of 2 Gaussians to the data with the below code (the plot.normal.components function is given at the end of my question):

if(!require("mixtools")) { install.packages("mixtools");  require("mixtools") }
snoq.k2 <- normalmixEM(snoq, k=2, maxit=100, epsilon=0.01)
plot(hist(snoq,breaks=101),col="grey",border="grey",freq=FALSE,
xlab="Precipitation (1/100 inch)",main="Precipitation in Snoqualmie Falls")
lines(density(snoq),lty=2)
sapply(1:2,plot.normal.components,mixture=snoq.k2)


The dotted line is the kernel density estimation (on the empirical pdf) and the two black curves are the fitted Gaussians. Now, I thought the distribution of snoq is similar to that of an exponential distribution, so my first instict here was to log transform the data and then investigate what happens if I try to fit a Mixture of Gaussians on log data rather than the raw data as the instructor did:

log_snoq <-  log(snoq)
log_snoq.k3 = normalmixEM(log_snoq, k=3, maxit=100, epsilon=0.01) # does not converge!
plot(hist(log_snoq,breaks=101),col="grey",border="grey",freq=FALSE,
xlab="Log Precipitation (1/100 inch)",main="Precipitation in Snoqualmie Falls")
lines(density(log_snoq),lty=2)
sapply(1:3,plot.normal.components,mixture=log_snoq.k3)


I chose k=3 in this case, just because the kernel density estimator showed 3 bumps on the dotted line (in the figure below), which I interpreted as three local maxima that can be modelled by 3 univariate Gaussians. However, during EM I get a warning that says WARNING! NOT CONVERGENT!, and the resulting Gaussian components look as below: (I do realise that the log transform does not produce a beautifully Gaussian-looking distribution, but the raw data itself looks less suitable to be modelled with a (mixture of) Gaussian(s) to me anyway.) My question is, in this particular example is it wrong to log-transform? Does this in general complicate the application of Mixture of Gaussians? Any recommendations / comments are welcome!

The function plot.normal.components:

plot.normal.components <- function(mixture,component.number,...) {
curve(mixture$lambda[component.number] * dnorm(x,mean=mixture$mu[component.number],
}

• I've added the 'plot.normal.components()' function. – ocram Feb 4 '15 at 8:32
• Sorry, it was already there, but at the bottom so probably not very visible – Zhubarb Feb 4 '15 at 8:33
• Oh, you are right! I let you decide which one to remove – ocram Feb 4 '15 at 8:36
• For a similar analysis as part of EDA for my dissertation research software, I wanted to automate mixture analysis as much as possible. Therefore, I've used mclust package to analytically determine number of mixture components and then pass to normalmixEM(), etc. You can take a look at that module. This worked well for my data, but when I've tried similar code for yours, it failed to detect the correct number of components (2 or 3) - it detected only one. Not sure why. – Aleksandr Blekh Feb 4 '15 at 15:11
• @AleksandrBlekh, did you try it on the raw or log(data)? Also, have you seen that in the original code on the link I gave, the instructor calculates the final loglikelihoods of different component numbers by 2-fold cross validation (fitting the model trained on half the dataset to the other half). – Zhubarb Feb 4 '15 at 15:29

• Unless I'm missing something, I don't fully understand why for some reason you've used different number of mixture components in the second and third blocks of code (3 vs. 2). However, looking at the output of summary() for the mixture object, it seems that the 2nd component is negligible, so I assume that you decided to ignore it, hence the change.
• Warning "NOT CONVERGENT!" can be alleviated in many cases by increasing the number of iterations (for example, I've used the value of 500 instead of the default 100 for my similar analysis). I've tried your code and it converges after 124 iterations, so maxit=200 should be enough.