r - ガウス分布とガンマ分布の混合

Question

ガウス分布とガンマ分布の混合物からコンポーネント分布パラメーターを見つけるために、R でいくつかのスクリプト/パッケージを探しています (Python もそうします)。これまで、R パッケージの "mixtools" を使用してデータをガウスの混合物としてモデル化してきましたが、ガンマとガウスを組み合わせることでより適切にモデル化できると思います。

ありがとう

Answer 1

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