What distribution can be closely (or precisely) fit to the "5 number summary" statistics?

Question

I'm programming a web tool (="I'm a stats ignoramous who drifted here from stackoverflow.com") that allows scientists to enter predictions about the 5-number-summary stats for a variable. Entry is done using the UI-metaphor of a box-plot.

I'd like to allow scientists to visualize their input as a PDF/CDF, but I have to select an underlying distribution.

• I'm looking for a distribution that is "as normal as possible" while still being able to fit well to the 5-number summary pinning down the ~1st, 25th, 50th, 75th and ~99th percentiles.
• I started with the 3-param skew-normal, but it obviously doesn't have enough DOF to perfectly (or even closely) fit to the 5-input parameters
• I'm interpreting 'min' and 'max' as 1st and 99th percentiles. I know this is sketchy, but the numbers entered are speculative predictions (="Don't worry, I'm not screwing up the interpretation of measured data")
• Simplicity is a virtue. Ideally the distribution would be easy to easy-ish to do numerical parameter-estimation with (closed form would be the nicest, ala http://www.johndcook.com/blog/2010/01/31/parameters-from-percentiles/ , but that's pretty much shooting for the moon, doing non-linear optimization or something is fine)
• I've started looking at distributions like GSN/CSN, etc from papers like http://www2.warwick.ac.uk/fac/sci/statistics/crism/research/2012/paper12-08/12-08w.pdf , but I'm not really sure I'm looking down the right family. Maybe skew-normal isn't the best place to start? I've also thought about things like Johnson distribution, which from the little I can find about it seems almost "designed to be fitted".

What distribution(s) should I be looking at?

Screenshot of the tool: failing to fit the median on a skew-normal distribution

Answer 1

user1448319's answer triggered the following thought in my brain. Do a natural cubic spline on the set of points of the form

$(x_p, \Phi^{-1}(p))$

where $x_p$ is the $100p$ percentile and $\Phi^{-1}(\cdot)$ is the quantile function of the normal distribution. Run the resulting interpolating spline function through the normal CDF and take the derivative to obtain the PDF. This procedure has the following properties:

• the resulting distribution matches the given percentiles exactly;
• the tails are normal;
• if the given percentiles actually match those of some normal distribution, the output is that normal distribution;
• the numerical computations are dead easy and give analytical expressions for the PDF;
• the generalization to other target distributions is obvious.

But the proof is in the pudding. Let me whip up some R code...

elicit_distribution <- function(x, p, qfun = qnorm, pfun = pnorm, dfun = dnorm, range_factor = 1, length.out = 1000, ...)
{
  fun <- splinefun(x, qfun(p), method = "natural", ...)
  cdfun <- function(x) pfun(fun(x, deriv = 0))
  from <- min(x) - range_factor*diff(range(x))
  to <- max(x) + range_factor*diff(range(x))
  xval <- seq(from, to, length.out = length.out)
  list(cdfun = cdfun
      ,pdfun = function(x) fun(x, deriv = 1)*dfun(fun(x, deriv = 0))
      ,quantfun = approxfun(cdfun(xval),xval)
      )
}

plot_elicited_distribution <- function(x, p, qfun = qnorm, pfun = pnorm, dfun = dnorm, range_factor = 0.1, lwd = 2, ylab = "PDF", ...)
{
  dist <- elicit_distribution(x,p,qfun,pfun,dfun)
  from <- min(x) - range_factor*diff(range(x))
  to <- max(x) + range_factor*diff(range(x))
  curve(dist$pdfun(x), from = from, to = to, lwd = lwd, ylab = ylab, ...)
  lineseg <- function(x,y,...)
    points(c(x,x),c(0,y),type = "l", lwd = lwd, ...)
  col <- function(i) c("red","green")[1+((i-1)%%2)]
  xval <- dist$quantfun(p)
  for(i in 1:length(xval))
  {
    points(x[i], dist$pdfun(x[i]), col = col(i), pch = 16)
lineseg(xval[i],dist$pdfun(xval[i]), col = col(i))
  }
}

x <- c(5, 15, 17, 25, 46)
p <- c(0.01, 0.25, 0.5, 0.75, 0.99)
plot_elicited_distribution(x,p)

(Solid points plotted on the PDF curve show given values. Lines show percentiles of the generated distribution.)

Aw, crap. Add one more property to the list:

• no guarantee of unimodality

Let's try a smoothing spline instead. Code as before, except in "elicit_distribution" replace

fun <- splinefun(x, qfun(p), method = "natural")

with

splineobj <- smooth.spline(x, qfun(p))
fun <- function(x, deriv) predict(splineobj, x, deriv)$y

That's a bit better. It's quite similar to the skew-normal plot you posted but it seems to have a different trade-off for awkward percentiles, resulting in a slightly better fit at the median and a slightly worse fit at the 25% point.

Answer 2

Why not just use something like a piecewise linear distribution?

Let's say a scientist gives you the values $x_{01}, x_{25}, x_{50}, x_{75}, x_{99}$ which correspond to the 1%, ..., 99% of the unknown underlying distribution. We want to make a distribution where there is 1% of the mass to the left of $y_{01}$, ..., and 99% of the mass to the left of $y_{99}$.

Let's call this distribution function $f$, i.e., $f(x_t) = y_t$.

Let's suppose that the distribution has a finite $x_{00}$ and $x_{100}$. Let's also suppose that we know what $x_{00}$ is. For now, let's pick something like $x_{00} = x_{01} - |x_{25} - x_{01}|$ (just so we have a specific value to do debugging with or something). I'll come back to this later.

Set $y_{00} = 0$. Set $y_{01}$ so that the area under the line segment from $(x_{00},y_{00})$ to $(x_{01},y_{01})$ is equal to 1% (i.e., so $\int_{x_{00}}^{x_{01}} f(x) = 0.1$). This gives you a value for $y_{01}$. Now find the value for $y_{25}$ so that the area under the line segment from $(x_{01},y_{01})$ to $(x_{25},y_{25})$ is equal to 25%-1%=24%. Do this again to find $y_{50}$, $y_{75}$, and $y_{99}$. No select the $x_{100}$ which gives you a total area of 100% under the piecewise linear function you've constructed. Now you have a distribution with exactly 1% of the mass to the left of the 1% value the expert told you, 25% to the left of the 25% value the expert told you, etc.

Now, look at your distribution. Pick a value of $x_{00}$ that makes sense. It might be clever to pick some measure that you want to minimize to give you an automatic selection of $x_{00}$. For example, you could minimize the total angle of your distribution (e.g., if your distribution is $f$, you could minimize $\int_{-\infty}^\infty {d^2 \over dx^2}f(x) dx$ which is just the sum of the angles of $f$ at each of $x_{00},...,x_{100}$).

This seems like the most naive approach to me, it is very flexible, and it has the added benefit of being non parametric so you don't have to estimate anything. I hope it's a good starting place.

Answer 3

You could achieve this based on Box-Cox transformation or other power transformation family (denpending on whether your random variable is strictly positive or not). First, your can assume the original unknown distribution is well-behaved (not from a mixed distribution). Then based on the Box-Cox transformation, the transformed distribution will be approximately normally distributed.

(1).Set the initial value of summary statistics for a normally distributed random variable. The initial values can be calculated by applying Box-Cox transformation to your reported summary statistics of the unknow distribution $X$. This will give you the initial values of $y_q$ and the initial transformation parameter $\lambda$.

(2). Simulate a normal random variable for the sample size of the study with the initial values from (1), therefore $y\sim Normal(\mu, \sigma^2)$. If you using quantiles in (1), then $\mu$ and $\sigma^2$ can be derived by using the formula of $\mu\pm v_q\sigma=y_q$, where $v_q$ is the theoretical quantile values for the normal distribution.

(3). Inverse the Box-Cox transformation $x=(y\lambda+1)^{1/\lambda}$ and calcualte summary statistics of sample mean, sample standard deviation or sample percentile ranges from inversed distribution of $x$.

(4). Minimize the sum of least-squares $\sum{\frac{\theta_i-O_i}{O_i}} $ to obtain the optimal estimates of normal random variable $Y$, where $\Theta$ is the vector of summary statistics from the inversed distribution, and $O$ is the vector of reported summary statistics from the unknown distribution.

(5). Substitute those optimal estimates into (2) and (3) to get the simulated distribution of the unknown.

(6). Go back to (2) and using different random seeds to simulate a new normal distribution.

I hope this helps.