Geek Answers Handbook

A mixture of gaussian and gamma distribution

I am looking for some script / package in R (Python will do as well) to find out the distribution parameters of components from a mixture of Gaussian and gamma distributions. I have so far used the mixools R package to simulate data as a mixture of Gaussians, but I think it can be better modeled by Gamma plus Gauss.

thanks

+4
source share
1 answer

Here is one possibility:

Define utility functions:

rnormgammamix <- function(n,shape,rate,mean,sd,prob) { ifelse(runif(n)<prob, rgamma(n,shape,rate), rnorm(n,mean,sd)) }

(This can be done a little more efficient ...)

 dnormgammamix <- function(x,shape,rate,mean,sd,prob,log=FALSE) { r <- prob*dgamma(x,shape,rate)+(1-prob)*dnorm(x,mean,sd) if (log) log(r) else r }

Generate fake data:

 set.seed(101) r <- rnormgammamix(1000,1.5,2,3,2,0.5) d <- data.frame(r)

Approach No. 1: bbmle package. Set the shape, speed, standard deviation on the logarithm scale, test logit scale.

 library("bbmle") m1 <- mle2(r~dnormgammamix(exp(logshape),exp(lograte),mean,exp(logsd), plogis(logitprob)), data=d, start=list(logshape=0,lograte=0,mean=0,logsd=0,logitprob=0)) cc <- coef(m1) png("normgam.png") par(bty="l",las=1) hist(r,breaks=100,col="gray",freq=FALSE) rvec <- seq(-2,8,length=101) pred <- with(as.list(cc), dnormgammamix(rvec,exp(logshape),exp(lograte),mean, exp(logsd),plogis(logitprob))) lines(rvec,pred,col=2,lwd=2) true <- dnormgammamix(rvec,1.5,2,3,2,0.5) lines(rvec,true,col=4,lwd=2) dev.off()

enter image description here

 tcc <- with(as.list(cc), c(shape=exp(logshape), rate=exp(lograte), mean=mean, sd=exp(logsd), prob=plogis(logitprob))) cbind(tcc,c(1.5,2,3,2,0.5))

The approach is reasonable, but the parameters are rather far - I think that this model is not very identifiable in this parameter mode (i.e. gamma and Gaussian components can be replaced)

 library("MASS") ff <- fitdistr(r,dnormgammamix, start=list(shape=1,rate=1,mean=0,sd=1,prob=0.5)) cbind(tcc,ff$estimate,c(1.5,2,3,2,0.5))

fitdistr gets the same result as mle2 , which means we're at a local minimum. If we start with true parameters, we will get something reasonable and close to true parameters.

 ff2 <- fitdistr(r,dnormgammamix, start=list(shape=1.5,rate=2,mean=3,sd=2,prob=0.5)) -logLik(ff2) ## 1725.994 -logLik(ff) ## 1755.458
+5
source

More articles:


All Articles