# Estimating a distribution based on three percentiles

What methods can I use to infer a distribution if I know only three percentiles?

For example, I know that in a certain data set, the fifth percentile is 8,135, the 50th percentile is 11,259, and the 95th percentile is 23,611. I want to be able to go from any other number to its percentile.

It's not my data, and those are all the statistics I have. It's clear that the distribution isn't normal. The only other information I have is that this data represents government per-capita funding for different school districts.

I know enough about statistics to know that this problem has no definite solution, but not enough to know how to go about finding good guesses.

Would a lognormal distribution be appropriate? What tools can I use to perform the regression (or do I need to do it myself)?

-
i've added the r tag so the R code is highlighted in my comment – mpiktas Jan 7 '11 at 13:50
For a detailed example of the same question (and its solution), see the duplicate thread at stats.stackexchange.com/questions/133129. – whuber Jan 12 '15 at 18:48

Using a purely statistical method to do this work will provide absolutely no additional information about the distribution of school spending: the result will merely reflect an arbitrary choice of algorithm.

You need more data.

This is easy to come by: use data from previous years, from comparable districts, whatever. For example, federal spending on 14866 school districts in 2008 is available from the Census site. It shows that across the country, total per-capita (enrolled) federal revenues were approximately lognormally distributed, but breaking it down by state shows substantial variation (e.g., log spending in Alaska has negative skew while log spending in Colorado has strong positive skew). Use those data to characterize the likely form of distribution and then fit your quantiles to that form.

If you're even close to the right distributional form, then you should be able to reproduce the quantiles accurately by fitting one or at most two parameters. The best technique for finding the fit will depend on what distributional form you use, but--far more importantly--it will depend on what you intend to use the results for. Do you need to estimate an average spending amount? Upper and lower limits on spending? Whatever it is, you want to adopt some measure of goodness of fit that will give you the best chance of making good decisions with your results. For example, if your interest is focused in the upper 10% of all spending, you will want to fit the 95th percentile accurately and you might care little about fitting the 5th percentile. No sophisticated fitting technique will make these considerations for you.

Of course no one can legitimately guarantee that this data-informed, decision-oriented method will perform any better (or any worse) than some statistical recipe, but--unlike a purely statistical approach--this method has a basis grounded in reality, with a focus on your needs, giving it some credibility and defense against criticism.

-
+1, great answer. – mpiktas Jan 6 '11 at 17:39
+1 You need more data and what you intend to use the results for deserve extra emphasis. – vqv Jan 6 '11 at 19:20
It sounds like there's lots of wisdom in your answer. I'll have to consult more with the people who posed me the problem about just what they want. Thank you for the links and the advice. – Mark Eichenlaub Jan 6 '11 at 19:29
@vqv Thanks. I made the changes you suggested. – whuber Jan 6 '11 at 19:33
@Mark Best of luck! – whuber Jan 6 '11 at 19:33

As @whuber pointed out, statistical methods do not exactly work here. You need to infer the distribution from other sources. When you know the distribution you have a non-linear equation solving exercise. Denote by $f$ the quantile function of your chosen probability distribution with parameter vector $\theta$. What you have is the following nonlinear system of equations:

\begin{align*} q_{0.05}&=f(0.05,\theta) \\ q_{0.5}&=f(0.5,\theta) \\ q_{0.95}&=f(0.95,\theta)\\ \end{align*}

where $q$ are your quantiles. You need to solve this system to find $\theta$. Now for practically for any 3-parameter distribution you will find values of parameters satisfying this equation. For 2-parameter and 1-parameter distributions this system is overdetermined, so there are no exact solutions. In this case you can search for a set of parameters which minimizes the discrepancy:

\begin{align*} (q_{0.05}-f(0.05,\theta))^2+ (q_{0.5}-f(0.5,\theta))^2 + (q_{0.95}-f(0.95,\theta))^2 \end{align*}

Here I chose the quadratic function, but you can chose whatever you want. According to @whuber comments you can assign weights, so that more important quantiles can be fitted more accurately.

For four and more parameters the system is underdetermined, so infinite number of solutions exists.

Here is some sample R code illustrating this approach. For purposes of demonstration I generate the quantiles from Singh-Maddala distribution from VGAM package. This distribution has 3 parameters and is used in income distribution modelling.

 q <- qsinmad(c(0.05,0.5,0.95),2,1,4)
points(p<-c(0.05, 0.5, 0.95), dsinmad(p, 2, 1, 4))


Now form the function which evaluates the non-linear system of equations:

 fn <- function(x,q) q-qsinmad(c(0.05, 0.5, 0.95), x[1], x[2], x[3])


Check whether true values satisfy the equation:

 > fn(c(2,1,4),q)
[1] 0 0 0


For solving the non-linear equation system I use function nleqslv from package nlqeslv.

 > sol <- nleqslv(c(2.4,1.5,4.3),fn,q=q)
> sol$x [1] 2.000000 1.000000 4.000001  As we see we get the exact solution. Now let us try to fit log-normal distribution to these quantiles. For this we will use the optim function.  > ofn <- function(x,q)sum(abs(q-qlnorm(c(0.05,0.5,0.95),x[1],x[2]))^2) > osol <- optim(c(1,1),ofn) > osol$par
[1] -0.905049  0.586334


Now plot the result

  plot(x,dlnorm(x,osol$par[1],osol$par[2]),type="l",col=2)


From this we immediately see that the quadratic function is not so good.

Hope this helps.

-
Great! Thanks for all the effort that went into this, mpiktas. I'm not familiar with R, but your code is explained well enough that I can still easily tell what you're doing. – Mark Eichenlaub Jan 7 '11 at 13:57
Thanks a lot for this example. I think there is 2 mistakes in ofn <- function(x,q) sum(abs(q-qlnorm(c(0.05,0.5,0.95),x[1],x[2]))^2). I propose ofn <- function(x) sum(abs(q-qlnorm(c(0.05,0.5,0.95),x[1],x[2],x[3]))^2) because q is not an input for ofn, and X[3] is missing. Regards – user25359 May 7 '13 at 23:00

For a lognormal the ratio of the 95th percentile to the median is the same as the ratio of the median to the 5th percentile. That's not even nearly true here so lognormal wouldn't be a good fit.

You have enough information to fit a distribution with three parameters, and you clearly need a skew distribution. For analytical simplicity, I'd suggest the shifted log-logistic distribution as its quantile function (i.e. the inverse of its cumulative distribution function) can be written in a reasonably simple closed form, so you should be able to get closed-form expressions for its three parameters in terms of your three quantiles with a bit of algebra (i'll leave that as an exercise!). This distribution is used in flood frequency analysis.

This isn't going to give you any indication of the uncertainty in the estimates of the other quantiles though. I don't know if you need that, but as a statistician I feel I should be able to provide it, so I'm not really satisfied with this answer. I certainly wouldn't use this method, or probably any method, to extrapolate (much) outside the range of the 5th to 95th percentiles.

-
Thanks for the advice. Re: lognormal - I could make the ratios of percentiles to median work out by subtracting 7077 from everything, then adding it back in at the end. How bad an idea would that be? – Mark Eichenlaub Jan 6 '11 at 9:04
Good point, that would give a 'shifted log-normal distribution'. The log-normal and the log-logistic are pretty similar in shape apart from the heavier tails of the latter, so you could try both and compare results. – onestop Jan 6 '11 at 11:40
Compare how? The shifted lognormal is guaranteed to fit the quantiles perfectly. Almost any three-parameter family will fit perfectly. How do you compare two perfect fits? – whuber Jan 6 '11 at 16:31
@whuber I meant compare the resulting predictions for the percentiles corresponding to other values – onestop Jan 6 '11 at 17:13
I'm missing something: what other values? The OP states that only three percentiles are available, nothing else. – whuber Jan 6 '11 at 17:38

Try the rriskDistributions package, and -- if you are sure about the lognormal distribution family -- use the command

get.lnorm.par(p=c(0.05,0.5,0.95),q=c(8.135,11.259,23.611))


which should solve your problem. Use fit.perc instead if you do not want to restrict to one known pdf.

-

About the only things you can infer from the data is that the distribution is nonsymmetric. You can't even tell whether those quantiles came from a fitted distribution or just the ecdf.

If they came from a fitted distribution, you could try all the distributions you can think of and see if any match. If not, there's not nearly enough information. You could interpolate a 2nd degree polynomial or a 3rd degree spline for the quantile function and use that, or come up with a theory as to the distribution family and match quantiles, but any inferences you would make with these methods would be deeply suspect.

-
Polynomials and splines are unlikely to be valid CDFs. – whuber Jan 6 '11 at 16:05
Good observation. In this case, the usual quadratic polynomial fails to work, but there are infinitely many quadratic splines to choose from (think Bézier) that should not have the same problem (though some might still require domain cropping). Similarly, it should be possible to find a suitable monotonic cubic spline. I am aware of spline algorithms that guarantee monotonicity, but am unable to find one just now, so I have to leave the matter at "pick something you like that works as cdf". – sesqu Jan 6 '11 at 17:46
You could go so far as to fit a monotonic spline (or whatever) to the logarithms of the quantiles, thereby obtaining something reasonable within the range of the quantiles. But this provides no help in fitting the tails beyond the two extreme quantiles. One should be reluctant to let such an important aspect of the fit be left to the accidental characteristics of the numerical fitting procedure. – whuber May 8 '13 at 13:17

The use of quantiles to estimate parameters of a priori distributions is discussed in the literature on human response time measurement as "quantile maximum probability estimation" (QMPE, though originally erroneously dubbed "quantile maximum likelihood estimation", QMLE), discussed at length by Heathcote and colleagues. You could fit a number of different a priori distributions (ex-Gaussian, shifted Lognormal, Wald, and Weibull) then compare the sum log likelihoods of the resulting best fits for each distribution to find the distribution flavor that seems to yield the best fit.

-
Any three-parameter distribution is guaranteed to fit three quantiles perfectly. Thus it makes sense to use this approach to fit only one or two parameters. It also does not make any sense to compare a one-parameter fit to a two-parameter fit (with a different family) based on likelihood alone. – whuber Jan 6 '11 at 16:04
@whuber, re: "Any three-parameter distribution is guaranteed to fit three quantiles perfectly". I hadn't realized that, so good to know! re: "It also does not make any sense to compare a one-parameter fit to a two-parameter fit (with a different family) based on likelihood alone." Ah yes, indeed; I failed to mention that one would have to apply some complexity correction (AIC, BIC, ...) if comparing fits to distribution flavors with different numbers of parameters. Thanks for pointing that out. – Mike Lawrence Jan 6 '11 at 21:10
I exaggerated a little bit, because I was thinking of two of the parameters being scale and location and the third comprising a wide range of shapes. Even so, most three-parameter families have sufficient flexibility to fit three percentiles provided they are all distinct. – whuber Jan 6 '11 at 23:11

You can use your percentile information to simulate the data in some way and use the R package "logspline" to estimate the distribution nonparametrically. Below is my function that employs a method like this.

calc.dist.from.median.and.range <- function(m, r)
{
## PURPOSE: Return a Log-Logspline Distribution given (m, r).
##          It may be necessary to call this function multiple times in order to get a satisfying distribution (from the plot).
## ----------------------------------------------------------------------
## ARGUMENT:
##   m: Median
##   r: Range (a vector of two numbers)
## ----------------------------------------------------------------------
## RETURN: A log-logspline distribution object.
## ----------------------------------------------------------------------
## AUTHOR: Feiming Chen,  Date: 10 Feb 2016, 10:35

if (m < r[1] || m > r[2] || r[1] > r[2]) stop("Misspecified Median and Range")

mu <- log10(m)
log.r <- log10(r)

## Simulate data that will have median of "mu" and range of "log.r"
## Distribution on the Left/Right: Simulate a Normal Distribution centered at "mu" and truncate the part above/below the "mu".
## May keep sample size intentionaly small so as to introduce uncertainty about the distribution.
d1 <- rnorm(n=200, mean=mu, sd=(mu - log.r[1])/3) # Assums 3*SD informs the bound
d2 <- d1[d1 < mu]                   # Simulated Data to the Left of "mu"
d3 <- rnorm(n=200, mean=mu, sd=(log.r[2] - mu)/3)
d4 <- d3[d3 > mu]                   # Simulated Data to the Right of "mu"
d5 <- c(d2, d4)                     # Combined Simulated Data for the unknown distribution

require(logspline)
ans <- logspline(x=d5)
plot(ans)
return(ans)
}
if (F) {                                # Unit Test
calc.dist.from.median.and.range(m=1e10, r=c(3.6e5, 3.1e12))
my.dist <- calc.dist.from.median.and.range(m=1e7, r=c(7e2, 3e11))
dlogspline(log10(c(7e2, 1e7, 3e11)), my.dist) # Density
plogspline(log10(c(7e2, 1e7, 3e11)), my.dist) # Probability
10^qlogspline(c(0.05, 0.5, 0.95), my.dist) # Quantiles
10^rlogspline(10, my.dist) # Random Sample
}

-