The ex-gauss distribution is well known for its ability to fit human reaction times distributions. I wanted to try something out, and was struck by the fact that it was so hard to find online expressions and source code. So I set out to do it myself, and defined the appropriate R functions.
The ex-gauss is, barebones, the sum of a normally (gaussian) distributed random variable with mean μ and variance σ, and an exponentially distributed random variable with rate parameter λ; assuming that the result is non-negative. So the simulation function is quite easy:
"rexgauss" <-
function(n,mu=0.5,sigma=.1,lambda=1)
rnorm(n,mu,sigma)+rexp(n,lambda)
The default values I chose are just to remind me about the range of sensible values for reaction time distributions. Next, the density function needs to be defined. Unfortunately, the only reference I had (Luce, 1986) contains an error, but it's not difficult to recalculate. The resulting function is
"dexgauss" <-
function (t, mu = 0.5, sigma = .1, lambda = 1,log=FALSE)
if(!log) lambda*exp(-t*lambda)*exp(mu*lambda+(lambda*sigma)^2/2) * pnorm(t,mu+lambda*sigma^2,sigma) else log(lambda)-t*lambda+mu*lambda+(lambda*sigma)^2/2 + pnorm(t,mu+lambda*sigma^2,sigma,log=TRUE)
Note that I included the option to have the log of the density returned, that's why it looks perhaps a bit messy.
Next, it would come in handy if we had the distribution function. This can actually be computed in terms of normal distributions (found by straightforward integration by parts of the density), so that we can use the efficient R routine pnorm. The resulting function is (not allowing for log transformed probabilities, maybe I should change that):
"pexgauss" <-
function(t, mu=0.5, sigma=0.1, lambda=1) pnorm(t,mu,sigma) - dexgauss(t,mu,sigma,lambda)/lambda
I could not find an analitical inverse of the distribution function, so for quantiles you'll have to create an approximation yourself (e.g., by using splinefun). Well, here's a function that implements the bisection method for this purpose, make sure you inspect the outcome:
qexgauss <-
function(p,mu=.5,sigma=.1,lambda=1){P=p;p=p[(not1 <- p<1) & (not0 <- p>0)]; f=function(t)pexgauss(t,mu,sigma,lambda)-p; b=a=p;v=sigma+1/lambda^2;a[]=mu-4*v;b=mu+10*v;fa=f(a);fb=f(b);c=co=b;co[]=0; while(max(abs(c-co)>1e-8)){co=c;c=(a+b)/2;fc=f(c); iac=fa*fc<0; a=ifelse(iac,a,c);fa=ifelse(iac,fa,fc); b=ifelse(iac,c,b);fb=ifelse(iac,fc,fb);};Q=P; Q[not1¬0]=(a+b)/2;Q[!not1]=Inf;Q[!not0]=-Inf; Q}
Because not all sets of parameter values are created equal—that is, not all combinations of parameter values correspond to sensible response time distributions—and the routines given here were designed on the premise that the results correspond to a sensible response time distributions, you have to test if the parameter values you choose result in a valid distribution. A simple test function does this for you (allowing for a small margin of error):
"valid.exgauss" <-
function (mu, sigma, lambda)
pexgauss(0, mu, sigma, lambda) > -1e-06
In most cases you'd also want to be able to fit the ex-gauss distribution to your data. To that end, you may use:
fit.exgauss.mom <- function(rt){lmb=1/(mean((rt-(m<-mean(rt)))^3)/2)^(1/3); c(mu = m-1/lmb, sigma = if((v<-var(rt))>1/lmb^2) sqrt(v-1/lmb^2) else 0, lambda = lmb)}
which returns method-of-moments estimates. A simple example of its use would be
rt = rexgauss(150); # simulate some response times
fit.exgauss.mom(rt); # get MoM estimates
Alternatively you may want to rely on asymptotic efficient estimators, and use
"fit.exgauss" <-
function (x, start = fit.exgauss.mom(rt), ...)
optim(start, function(theta, y) -sum(dexgauss(y, mu = theta[1],
sigma = theta[2], lambda = theta[3], log = TRUE)), y = x, ...)
which returns maximum likelihood estimates. An example use is
fit.exgauss(rexgauss(150))
This routine uses the Nelder-Mead simplex optimization algorithm as the default, which is quite stable (it evens seems to work reasonable well with only 50 observations). The downside is that you don't get a Hessian, or standard errors for that matter. Instead of Nelder-Mead you can use any of the other optimization algorithms that optim provides, by specifying extra arguments which are directly passed to the optim routine (hence, see the optim help file for all available arguments). For instance, you can use the BFGS algorithm with lower bounds to avoid warning messages about NaN's and non-finite objective function values:
fit.exgauss(rexgauss(150), method="L-BFGS", lower = 0.00001, hessian=TRUE) # also return a hessian
Here, lower is set to a small number to indicate all parameters should be greater than 0. It is not set equal to zero because if we would have, optim would actually try to set sigma to zero which results in an error because the ex-gaussian density is not defined for this parameter value. The argument hessian is set to TRUE to have the Hessian returned after convergence. The inverse of the Hessian is an estimate of the covariance matrix of the parameter estimates from which you can calculate the parameter estimate standard errors.
Alternatively, you might want to try nlm instead of optim. To make use of the full power of the gradient based optimization in nlm, derivatives are needed. So, I created (...well, R & I did it together... we both couldn't have done it without each other) a negative log-likelihood function that also returns the gradient in the way that nlm likes it: First the log of dexgauss with gradients (this is R's part of the work)
"ddexgauss.log" <-
deriv(~log(lambda*exp(-t*lambda)*exp(mu*lambda+(lambda*sigma)^2/2) * pnorm((t-mu-lambda*sigma^2)/sigma)),
c('mu','sigma','lambda'), function(t, mu=0.5, sigma=0.1, lambda=1){})
Next the negative log-likelihood function (this is where I came in):
"nloglike.exgauss" <-
function (theta, y, gradients = FALSE)
{
.tmp = ddexgauss.log(y, theta[1], theta[2], theta[3]);
.val = -sum(.tmp);
if (gradients)
attr(.val, "gradient") <- colSums(-attr(.tmp, "gradient"));
.val;
}
Note that it doesn't return the gradients by default, so in nlm you will need to add the parameter 'gradients=TRUE'. As an example of its use:
nlm(nloglike.exgauss,c(.5,.1,1), hessian=TRUE, y=rexgauss(500), gradients=TRUE)
Notice that I simulated 500 responses in this example. This is because nlm is far more sensitive to strange behavior of the likelhood function (e.g., due to invalid parameter values) than the Nelder-Mead simplex optimizer, and tends to exit more easily with an error message issued, than the Nelder-Mead algorithm. This can mostly be alleviated by specifying the maximum step size in nlm (e.g., stepmax=0.1, please read help(nlm)). The use of fit.exgauss may therefore be more convenient. However, in addition to returning the Hessian matrix upon request, nlm does use far fewer function evaluations.
That's all folks!
Update: See also this post for a routine implementing the QMLE method of Heathcote et al. (2002).
Update: Please read this report on a comparison between the estimates obtained with these routines and other software.
Labels: density, distributions, ex-gauss, fitting, maximum likelihood, simulation
