Using R to estimate probability distribution parameters using quantiles or HPD intervals

This post is written by Jordan Douglas. If you know some of the quantiles of a probability distribution (eg. `P(X < 5) = 0.05`, `P(X < 10) = 0.95`) and you know the name of the distribution but you do not know its parameters, then you can use R to estimate them. This can be done using the cumulative distribution functions (pnorm for normal, plnorm for lognormal, pexp for exponential etc.). Parameter values which yield a small difference between observed and calculated cumulative probabilities can be found using the optim function. For example, suppose that you want to find mu and sigma parameters of a lognormal distribution such that `P(X < 5) = 0.05` and `P(X < 10) = 0.95`. We start at `mu = 1` and `sigma = 1`. These initial values are a poor fit. Next we generate an error function, in this case we minimize the sum of squared differences between calculated and observed cumulative probabilities. Finally, we call optim to find estimates for sigma and mu: Thus, the estimates for mu and sigma are 1.956 and 0.211 respectively, which have a sum of squared residuals of 5.5e-11.