聞いたこともないけど,要素数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)

*追記(1/28)
標本サイズ,母数,初期値次第では結構うまくいく場合もあるようだ。たとえば,
set.seed(1)
dat<-c(rexgauss(500,1000,500,100),rexgauss(500,10000,500,100))
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)