r - Mixture of Gaussian and Gamma distribution -

i'm looking script/package in r (python too) find out component distribution parameters mixture of gaussian , gamma distributions. i've far used r package "mixtools" model data mixture of gaussians, think can better modeled gamma plus gaussian.

thanks

here's 1 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 made little bit 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 #1: bbmle package. fit shape, rate, standard deviation on log scale, prob on 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() 

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 fit reasonable, parameters far off -- think model isn't identifiable in parameter regime (i.e., gamma , gaussian components can swapped)

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 same result mle2, suggests we're in local minimum. if start true parameters reasonable , near 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