# Can anyone explain quantile maximum probability estimation (QMPE)?

Heathcote, Brown & Mewhort (2002, PDF) present Quantile Maximum Probability Estimation (originally termed Quantile Maximum Likelihood Estimation but later corrected) as a method of fitting distributional data, and find that it outperforms the more traditional Continuous Maximum Likelihood Estimation approach, at least in the case of fitting data to the ex-Gaussian (though see also this paper, which shows that this benefit generalizes to other distributions as well).

I'm trying to understand the actual steps to achieving QMPE. I understand that one first specifies increasing and equidistant quantile probabilities, then uses these to obtain the quantile values (q) in the observed data corresponding to these probabilities. I also understand that these observed quantile values are then used in counting the number of observations between each quantile (N). But this is where I'm stuck. Presumably one searches through the parameter space of whatever a priori model one assumes generated the data, searching for a parameter set that maximizes the joint probability of q and N. However, I don't know how, given a set of candidate parameters, this joint probability is computed.

Not having a strong math background, I think much better in code, so if someone could help me figure out what comes next, I'd be greatly appreciative. Here's the beginning of an attempt to fit some data to an ex-Gaussian:

    #generate some data to fit
true_mu = 300
true_sigma = 50
true_tau = 100
my_data = rnorm(100,true_mu,true sigma)+rexp(100,1/true_tau)

#select some quantile probabilities; estimate quantiles and inter-quantile counts from the observed data
quantile_probs = seq(.1,.9,.1) #or does it have to be seq(0,1,.1) ?
q = quantile( my_data, probs = quantile_probs, type = 5 ) #Heathcote et al apparently use type=5 given their example
N = rep( NA , length(q)-1 )
for( i in 1:( length(q)-1 ) ){
N = length( my_data[ (my_data>q[i]) & (my_data<=q[i+1]) ] )
}

#specify some candidate parameter values to assess (normally done as part of an iterative search  using an optimizer like optim)
candidate_mu = 350
candidate_sigma = 25
candidate_tau = 30
#given these candidates, what next?

-