QML Ex-gauss estimation

A while ago, I wrote about the Ex-Gauss distribution and posted some (rough) R-code for fitting the Ex-Gauss to data, using maximum likelihood. Now the Ex-Gauss can be quite hard to fit, especially with few observations (just try fit.exgauss(rexgauss(20)) and see the number of warnings). To mitigate this problem, Andrew Heathcote et al. (2002), proposed to use a modified maximum likelihood procedure, what they call "quantile maximum likelihood". Judith Dirk asked me if I had R code for that procedure as well. I didn't, so I wrote some. The code below implements this procedure for the Ex-Gauss distribution. You need to hand the routine two vectors: The vector of observations rt, and a vector p that specifies the percentiles that you want to use for quantile estimation. Although, I'm not sure about the statistical properties of this estimator, there seems to be some justification for asymptotic ML (Heathcote & Browne (2004) Psychon. Bul. and Rev.)

# This function fits the Ex-Gauss distribution using the Quantile
# Maximum Likelihood method of Heathcote (2002). It relies on
# the R routines defined on
# /rpages/2007/07/ex-gaussian-distribution-for-reaction.html
#
# example use:
# rt = rexgauss(150);
# p = seq(0, 1, len=13) # should contain 0 and 1!
# fit.exgauss.QMLE(rt, p)
#
function(rt, p, start =fit.exgauss.mom(rt), ...)
{
p = sort(unique(c(0,1,p)));
RT = sort(rt);
N = diff(p)*length(rt);
I = p[-c(1,length(p))]*length(rt)+ 0.5;
Im = trunc(I);
Ip=trunc(Im+1);
qhat = c(-Inf,RT[Im] + (RT[Ip]-RT[Im])*(I-Im), Inf);
f=function(theta, q, N) {
P=pexgauss(qhat,mu=theta[1],sigma=theta[2],lambda=theta[3]);
P[1]=0;
-sum(N*log(diff(P)));
}
optim(start, f, q=qhat, N=N, ...);
}

Update (Sep 18, 2008)

Because of stability issues it is advisable to restart the estimation procedure a number of times and retain the fit with the smallest objective function value. Furthermore, sometimes a different parameterization may be desirable. In particular, often θ = 1/λ is used for the exponential part of the ex-gaussian distribution. To facilitate the use of the above routine, and the earlier routines (including those for method-of-moments and ML), I have defined a new function, fitExgauss that will automate repeated fitting the ex-gauss using either standard ML or QMLE, and returns both parameterizations plus standard errors. The routine is defined in the file exgauss.R that you can download here: exgauss.R

The use is as follows:

source("exgauss.r");
rt = rexgauss(150, mu=500, sigma=100, lambda=1/80); # simulate response times
fitExgauss(rt); # gives ML estimates
fitExgauss(rt, seq(0,1, len=9)); # gives QMLE estimates using 9 bins

You can specify the number of repeated fitting rounds (nrepeat argument), different starting values (set start parameter to either a 3-vector of starting values or a function that accepts the array of response times and returns a 3-vector of starting values), or any of the optional arguments of optim (see ?optim for more).

Update (September, 2009): Update: Please read this report on a comparison between the estimates obtained with these routines and other software.

uploaded_files/exgauss.R

EZ2: An extension of the EZ-diffusion model for Response Time and Accuracy

Yesterday, I submitted a paper to the Journal of Mathematical Psychology. In the paper I promised to have an R package available for download from the internet. So I'd better live up to that promise, otherwise the reviewers might nog accept the paper! Here it is!. (A version that is installable on Mac OS X with the package installer option 'Local Source Package' is here.) As an alternative to ex-Gauss modeling of response times, this code fits diffusions that model the way information accumulation takes place in the brain (or at least, that's what we think).

For Windows users, you may need to install it with the rcmd command line utility, but you would have to have Perl installed. I'm trying to make it available as a regular zip file soon. (added: here it is.) Meanwhile you might simply try this source file which comes without documentation.

Ahum... the package itself is hardly documented at this point. The documentation essentially boils down to the following:

## create some data (theoretical values, not simulated)
# clearly, you would like to fit real data!
A = seq(.08,.13,len=6)
X2 = data.frame(A=A)
X2$vrt0 = sapply(A, function(a) EZ2.vrt(.1,.05,a))
X2$pe0 = sapply(A, function(a) EZ2.pe(.1,.05,a))
X2$vrt1 = sapply(A, function(a) EZ2.vrt(.2,a-.05,a))
X2$pe1 = sapply(A, function(a) EZ2.pe(.2,a-.05,a))

X2 = as.data.frame(X2)

# now pretend that X2 is the data frame that you have computed from real data

## fit an EZ2 model on each row
# method 1:
EZ2batch(c(v0=.11,v1=.21,z=.05,a=.09),
vrt0 ~ EZ2.vrt(v0,z,a),
pe0 ~ EZ2.pe(v0,z,a),
vrt1 ~ EZ2.vrt(v1,a-z,a),
pe1 ~ EZ2.pe(v1, a-z, a), data=X2)

# method 2 (eventually less typing):
mdl <- list( vrt0 ~ EZ2.vrt(v0,z,a),
pe0 ~ EZ2.pe(v0,z,a),
vrt1 ~ EZ2.vrt(v1,a-z,a),
pe1 ~ EZ2.pe(v1, a-z, a)
)
EZ2batch(c(v0=.11,v1=.21,z=.05,a=.09), mdl, data=X2)

You may also want to try the webapplication. Way more difficult I think, but with a graphical user interface.

EZ2_1.0.tar.gz
EZ2_1.0_R_i386-apple-darwin8.10.1.tar.gz
EZ2_1.0.zip