Смесь гауссовского и гамма-распределения

Question

Я ищу некоторый скрипт/пакет в R (Python тоже сделает), чтобы узнать параметры распределения компонентов из смеси гауссовских и гамма-распределений. Я до сих пор использовал пакет R «mixools», чтобы моделировать данные как смесь гауссианцев, но я думаю, что это может быть лучше смоделировано Gamma plus Gaussian.

Благодаря

Answer 1

Вот одна возможность:

Определение функции полезности:

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

(Это может быть немного более эффективным ...)

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 
} 

Генерация поддельные данные:

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

Подход №1: bbmle упаковка. Установите форму, скорость, стандартное отклонение по шкале журнала, пробную шкалу логита.

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

Подгонка является разумным, но параметры довольно далеко - я думаю, что эта модель не очень сильно идентифицируемых в этом параметре режима (т.е. гамма и гауссовой компоненты могут быть заменены)

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 получает тот же результат, как и mle2, который предполагает, что мы в локальном минимуме. Если мы начнем с истинных параметров, получим на что-то разумное и близкое к истинным параметрам.

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