# Is this a reasonable approach to fitting distributions?

Take the task of fitting an a priori distribution like the ex-Gaussian to a collection of observed human response times (RT). One method is to compute the sum log likelihood of each observed RT given a set of candidate ex-Gaussian parameters, then try to find the set of parameters that maximizes this sum log likelihood. I wonder if this alternative approach might also be reasonable:

1. Select a set of equidistant quantile probabilities, e.g.:

qps = seq( .1 , .9 , .1 )

2. For a given set of candidate ex-Gaussian parameters, estimate the quantile RT values corresponding to qps, e.g.:

sim_dat = rnorm( 1e5 , mu , sigma ) + rexp( 1e5 , 1/tau )
qrt = quantile( sim_dat , prob = qps )

3. For each sequential interval between the thus-generated quantile RT values, count the number of observations falling into that interval, e.g.:

obs_counts = rep( NA , length(qrt)-1 )
for( i in 1:(length(qrt)-1) ){
obs_counts[i] = length( obs_rt[ (obs_rt>qrt[i]) & (obs_rt<=qrt[i+1]) ] )
}

4. Compare these observed counts to the expected counts:

exp_counts = diff(range(qps)) * diff(qps)[1] * length(obs_rt)
chi_sq = sum( (( obs_counts - exp_counts )^2 )/exp_counts )

5. Repeat steps 2-4, searching for candidate parameter values that minimize chi_sq.

Is this approach a reasonable alternative to the more standard maximum likelihood estimation procedure? Does this approach already have a name?

Note that I use the example of an ex-Gaussian purely for illustrative purposes; in practice I'm playing with using the above approach in a rather more complicated context (e.g. fitting data to a stochastic model that yields multiple distributions, each with a different number of expected observation count). The purpose of this question is to ascertain whether I've re-invented the wheel as well as if anyone can pick out any problematic features of the method.

-

One problematic feature is that there may be a continuum of optimal solutions. In most settings the quantiles are continuous functions of the parameters. When the distributions are continuous, almost surely there will be positive intervals between the data values. Suppose your objective function is optimized by a particular parameter value whose quantiles do not coincide exactly with any of the data: that is, they lie in the interiors of the intervals determined by the nearby data values. (This is an extremely likely event.) Then small changes in the parameter value will move the quantiles slightly, to remain within the same intervals, thereby leaving the chi-squared value unchanged because none of the counts changes. Thus the procedure doesn't even pick out a definite set of parameter values!

Another problematic feature is that this procedure apparently provides no way to obtain estimation errors for the parameters.

Another problem is that you do not know even the most basic properties of this estimator, such as its amount of bias.

-

What you are proposing is called quantile matching, though the way you propose to do it will be exhausting. The ex-Gaussian distribution can be found in the package gamlss.dist with quantiles as qexGAUS etc.; it uses nu where you use tau.

A similar quantile matching method can be used in the function fitdist in the package fitdistrplus using method="qme". The package is mentioned in the answer linked by bill_080. One difference is that it only matches as many quantiles as there are parameters (three in this case).

The following seems to work more or less: it simulates some data points from a particular ex-Gaussian distribution and then tries to estimate the parameters using quantile matching, and then draws some graphs. It needs a rough estimate of the parameters to work.

library(fitdistrplus)
library(gamlss.dist)

set.seed(1)
sim_size <- 1000
Gm <- 10 # mean of Gaussian
Gs <- 2  # sd of Gaussian
Em <- 5  # mean of exponential
sim_dat <- rnorm( sim_size , Gm , Gs ) + rexp( sim_size , 1/Em )

fit_qme <- fitdist(sim_dat, "exGAUS", method="qme",
start=c(mu=15, sigma=1, nu=3),
probs=c(0.2,0.5,0.8)               )
fit_qme
plot(fit_qme)


In this example and with this seed, the estimates are

> fit_qme
Fitting of the distribution ' exGAUS ' by matching quantiles
Parameters:
estimate
mu    9.859207
sigma 1.753703
nu    5.049785


By comparison a maximum likelihood estimate method using the same function might look something like

fit_mle <- fitdist(sim_dat, "exGAUS", method="mle",
start=c(mu=15, sigma=1, nu=3)      )


and produce something like

> fit_mle
Fitting of the distribution ' exGAUS ' by maximum likelihood
Parameters:
estimate Std. Error
mu    9.938870  0.1656315
sigma 2.034017  0.1253632
nu    5.007996  0.2199171

-
I checked out the source for the qmedist function in fitdistrplus and it appears that that algorithm compares the observed versus expected quantile values given a set of quantile probabilities. What I've implemented above compares the observed versus expected density of data between a given set of quantiles. – Mike Lawrence May 23 '11 at 22:51
@Mike: You may be right, but I would still describe both as a form "quantile matching", albeit with different types of loss functions. – Henry May 23 '11 at 23:02
@Mike @Henry This subtle difference, though, is crucial: binning by sample quantiles has legitimate theoretical justification whereas dynamic re-binning by quantiles determined by the parameters does not and suffers from several problematic behaviors. – whuber May 24 '11 at 16:32

Take a look at the QQ-Plot (under my answer) in the following link:

Need help identifying a distribution by its histogram

-
A QQ-plot compares observed and expected quantile values, whereas what I've implemented above compares the observed versus expected density of data between a given set of expected quantiles. I edited my original question to clarify that I'm not seeking alternative fitting approaches, I'm more seeking input whether there are any obvious pitfalls to the one I've implemented, and also whether it is an already named method. – Mike Lawrence May 23 '11 at 23:01