In reply to Umberto.

Thanks for the comment. I passed it on to a colleague who knows about Bayesian analysis and he sent me the following comments and code regarding solving this problem with PROC MCMC in SAS:

Least square and MLE in normal likelihood are the same thing, so we can use a normal likelihood function, say, in NLMIXED:

proc nlmixed data=SampleQntls;
parms alpha 2 beta 2 s2 0.001;
mu = cdf("Gamma", q, alpha, beta);
model p ~ n(mu, s2);
run;

which gives the same estimates as alpha and beta, 2.7164 and 1.700.

Now we can do the same thing in PROC MCMC, and use noninformative priors, flat on alpha and beta, Jeffreys’ on the variance. (Well, I don’t know if uniform on alpha and beta are noninformative in this case because they enter the model via the CDF function, but the Jeffreys’ prior on variance in a normal model is more of a safe bet).

proc mcmc data=SampleQntls nmc=50000 seed=17073;
parms alpha 2 beta 2;
parms s2 0.001 / slice;
prior alpha beta ~ uniform(0, 10);
prior s2 ~ general(-log(s2));
mu = cdf("Gamma", q, alpha, beta);
model p ~ n(mu, var=s2);
run;

This gives reasonable estimates to alpha and beta:

Standard
Parameter N Mean Deviation 95% HPD Interval

alpha 50000 2.7237 0.3567 1.9922 3.4810
beta 50000 1.7331 0.2785 1.1881 2.2261
s2 50000 0.00112 0.00203 0.000062 0.00360

If one is concerned with mixing due to the Jeffreys’ prior, we can use an inverse-gamma prior that mimics the Jeffreys’, with mean pushed very close to zero:

proc mcmc data=SampleQntls nmc=50000 seed=17073;
parms alpha 2 beta 2;
parms s2 0.001;
prior alpha beta ~ uniform(0, 10);
prior s2 ~ igamma(shape=2, scale=0.0001);
mu = cdf("Gamma", q, alpha, beta);
model p ~ n(mu, var=s2);
run;

which effectively gave similar point estimates.

Standard
Parameter N Mean Deviation 95% HPD Interval

alpha 50000 2.7226 0.2028 2.3291 3.1454
beta 50000 1.7079 0.1393 1.4402 2.0068
s2 50000 0.000316 0.000247 0.000071 0.000737