読者です 読者をやめる 読者になる 読者になる

草薙の研究ログ

英語教育関係。でも最近は統計(特にR)ネタが中心。

混合正規指数合成分布モデル(?)を最尤法で…

r 統計

聞いたこともないけど,要素数2の混合正規指数合成分布(ex-Gaussian)モデルというのを考えてみる。ま,2つの異なる認知プロセスが混合したときの反応時間の分布だとか,そんなそれっぽいことを考えてみる。そんなものは多分ない。

ま,でもこの確率密度関数は,

dmixexgauss<-function(x,lambda,mu1,sigma1,tau1,mu2,sigma2,tau2){
	y<-lambda*(1/tau1)*exp(mu1/tau1+(sigma1^2)/(2*tau1^2)-x/tau1)* pnorm(x,mu1+(1/tau1)*sigma1^2, sigma1)
+(1-lambda)*(1/tau2)*exp(mu2/tau2+(sigma2^2)/(2*tau2^2)-x/tau2)* pnorm(x,mu2+(1/tau2)*sigma2^2,sigma2)
return(y)
	}

こうなはず。これで尤度関数を作って,bbmleパッケージのmle2とかで最尤推定してみれば…と思ったんだけど,やっぱ上手く推定できなかった。methodもBFGSとかいろいろ試したけど,ううむ。なんと7母数の確率密度関数だもんな。俺の間違いってのもあるし。

ま,でもこの自作の確率密度関数自体はあってるぽい。暇ができたらちょっとやってみよう。

library(retimes)
set.seed(1)
dat<-c(rexgauss(400,1000,300,800),rexgauss(100,4000,800,500))
hist(dat,main="",breaks=20,col="lightblue",freq=F,ylim=c(0,.0006),xlim=c(0,8000))
x<-seq(0,8000,1)
lines(x,dmixexgauss(x,.80,1000,300,800,4000,800,500),lty=2,col=2)


f:id:kusanagik:20170127211001p:plain


*追記(1/28)
標本サイズ,母数,初期値次第では結構うまくいく場合もあるようだ。たとえば,

#あからさまな例
set.seed(1)
dat<-c(rexgauss(500,1000,500,100),rexgauss(500,10000,500,100))

#mle2で最尤推定
#対数尤度を返す関数
lf<-function(l,m1,s1,t1,m2,s2,t2){
-sum(log(dmixexgauss(dat,l,m1,s1,t1,m2,s2,t2)))
}

#適当な初期値で最尤推定
fit<-mle2(lf,start=list(l=.01,m1=1000,s1=100,t1=100,m2=10000,s2=100,t2=100))
summary(fit)

Maximum likelihood estimation

Call:
mle2(minuslogl = lf, start = list(l = 0.01, m1 = 1000, s1 = 100, 
    t1 = 100, m2 = 10000, s2 = 100, t2 = 100))

Coefficients:
     Estimate Std. Error  z value     Pr(z)    
l  5.0001e-01 1.5811e-02  31.6234 < 2.2e-16 ***
m1 8.8352e+02 5.3678e+01  16.4598 < 2.2e-16 ***
s1 4.2277e+02 2.8256e+01  14.9618 < 2.2e-16 ***
t1 2.5258e+02 5.1428e+01   4.9114 9.043e-07 ***
m2 9.8383e+03 5.6083e+01 175.4242 < 2.2e-16 ***
s2 4.4467e+02 2.8829e+01  15.4245 < 2.2e-16 ***
t2 2.4854e+02 5.3408e+01   4.6536 3.263e-06 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

-2 log L: 16646.24 

coef(fit)
vcov(fit)
AIC(fit)