3
$\begingroup$

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. enter image description here

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

enter image description here

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],
                sd=mixture$sigma[component.number]), add=TRUE, ...)
}
$\endgroup$
  • $\begingroup$ I've added the 'plot.normal.components()' function. $\endgroup$ – ocram Feb 4 '15 at 8:32
  • $\begingroup$ Sorry, it was already there, but at the bottom so probably not very visible $\endgroup$ – Zhubarb Feb 4 '15 at 8:33
  • $\begingroup$ Oh, you are right! I let you decide which one to remove $\endgroup$ – ocram Feb 4 '15 at 8:36
  • $\begingroup$ 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. $\endgroup$ – Aleksandr Blekh Feb 4 '15 at 15:11
  • $\begingroup$ @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). $\endgroup$ – Zhubarb Feb 4 '15 at 15:29
1
$\begingroup$

Some comments, advice and answers:

  • 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.
  • I don't think that it's wrong (harmful) to log-transform in this case and it seems that it is even beneficial, as it exposes the structure of the data. However, we definitely need to be careful with log transformations: https://stats.stackexchange.com/a/130275/31372.
$\endgroup$
  • $\begingroup$ Thanks Aleksandr, just to reconfirm then, you agree: 1) log transforming the data (before running mixt of Gaussians) is a better (or at least justifiable) design choice given the e-pdf of this dataset; 2) there is no side-effect in terms of EM convergence when we operate on log data. (and the fact that I get a non-convergence warning is kind of a red herring) $\endgroup$ – Zhubarb Feb 4 '15 at 12:47
  • 1
    $\begingroup$ @Zhubarb: You're welcome. Well, I can confirm #1 (log transformation) for this particular case (data set) and to the best of my knowledge. However, I cannot confirm #2 (the absence of side effects). It makes sense IMHO to take a look at the E-M algorithm, just to be sure. Or to wait for other people's feedback. $\endgroup$ – Aleksandr Blekh Feb 4 '15 at 12:55
  • 1
    $\begingroup$ Re #2: the first thing I thought was potential loss of numerical precision due to log transform and that causing problems in EM convergence - but then again log is a monotonic function and that should actually not be an issue. Practically, it should not be different than finding MLE using loglikelihood. $\endgroup$ – Zhubarb Feb 4 '15 at 12:57
  • 1
    $\begingroup$ @Zhubarb: I believe you :-). I was curious enough to perform a brief research and it seems that you're right (for example, see this paper and this dissertation). Some researchers (see this paper) even introduce log transformation into EM algorithm in order to improve certain aspects (in this case, for better classification). $\endgroup$ – Aleksandr Blekh Feb 4 '15 at 13:24
  • $\begingroup$ That is cool, I will check the references out, thank you :) . $\endgroup$ – Zhubarb Feb 4 '15 at 13:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.