Estimating a distribution based on three percentiles

Question

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)?

Answer 1

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.

Answer 2

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)
 plot(x<-seq(0,2,by=0.01), dsinmad(x, 2, 1, 4),type="l")
 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 the function nleqslv from package nleqslv.

 > 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)
  lines(x,dsinmad(x,2,1,4))
  points(p,dsinmad(p,2,1,4))

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

Hope this helps.

Answer 3

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.

Answer 4

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.

Answer 5

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.

Answer 6

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.

Answer 7

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 
}