# Seeking a distribution, perhaps uncommon, consistent with two data points and expert constraints?

I am trying to state a prior distribution for a Bayesian meta-analysis.

I have the following information about a random variable:

1. Two observations: 3.0, 3.6
2. a scientist who studies the variable has told me that $P(X<2)=P(X>8)=0$, and that values as high as 6 have nonzero probability.

I have used the following approach to optimization (the mode of log-N = $e^{\mu-\sigma^2)}$:

prior <- function(parms, x, alpha) {
a <- abs(plnorm(x[1], parms[1], parms[2]) - (alpha/2))
b <- abs(plnorm(x[2], parms[1], parms[2]) - (1-alpha/2))
mode <- exp(parms[1] - parms[2]^2)
c <- abs(mode-3.3)
return(a + b + c)
}
v = nlm(prior,c(log(3.3),0.14),alpha=0.05,x=c(2.5,7.5))
x <- seq(1,10,0.1)
plot(x, dlnorm(x, v$estimate[1], v$estimate[2]))
abline(v=c(2.5,7.5), lty=2) #95%CI


In the figure, you can see the distribution that this returns, but I would like to find something more like the red lines I have drawn in.

This provides the same shaped distribution using the lognormal, gamma, or the normal, and it results in a distribution with $P(X=5)<0.05$ and $P(X=6)< 0.01$, i.e.:

 plnorm(c(5,6), v$estimate[1],v$estimate[2])


Can anyone suggest alternatives? I would prefer to stick with a single distribution rather than a mixture.

Thanks!

-
In your question you have $P(X = 2) = 0$ and $P(X = 8) = 0$. Do you mean: $P(X < 2) = 0$ and $P(X > 8) = 0$?? –  M. Tibbits Dec 12 '10 at 18:32
@David After reading this question over many, many times I still don't understand it. Are you using the two observations and the scientist's opinion to estimate a prior for a Bayesian analysis? Is your prior going to be based on the scientist's opinion only and then you want to update it with the observations? What is the distinction between "unsurprising" and having zero probability? –  whuber Dec 15 '10 at 17:24
@David I appreciate your care in eliciting the information. But assigning specific values of 0.05 and 0.01 to those probabilities is questionable. That's not your fault; it's just how things are. We can't expect people to pin down probabilities that well (that's what data are for). Maybe you would like to represent those with hyperpriors :-)? –  whuber Dec 15 '10 at 17:54
@whuber Thank you for following up. Now I understand the concept of the maximum entropy solution that you previously suggested, although it would require some learning on my end before I would want to apply it. Your alternative, maximizing variance, sounds sufficiently consistent with my objective and straightforward for me to implement. Thanks again. –  David Dec 15 '10 at 19:37
@David Be aware that a problem with maximizing variance is that the solution will be a discrete distribution. I expect it to concentrate 0.01 probability at X=8, 0.04 at X=6, 157/300 at X=2, and the rest (32/75) at X=5. (This variance equals 2.59.) –  whuber Dec 15 '10 at 20:38

If, given an answer to my comment above, you wish to bound the range of the distribution, why not simply fit a Beta distribution where you rescale to the unit interval? In other words, if you know that the parameter of interest should fall between $[2 , 8]$, then why not define $Y = \frac{X - 5}{6} + \frac{1}{2} = \frac{X - 2}{6}$. Where I've first centered the interval on zero, divided by the width so that Y will have a range of 1, and then added $\frac{1}{2}$ back so that the range of Y is $[0,1]$. (You can think of it either way: directly from $[2,8] \rightarrow [0,1]$ or from $[2,8] \rightarrow [-\frac{1}{2},\frac{1}{2}] \rightarrow [0,1]$, but I thought the latter might be easier at first).

Then, with two data points, you could fit a beta posterior with a uniform beta prior?

-

What about the Kumaraswamy distribution, which has the following pdf:

$$f(x; a,b) = a b x^{a-1}{ (1-x^a)}^{b-1}$$ for $a>0$, $b>0$, $0 < x < 1$. This distribution can be rescaled to have the required support.

-

Since the log-normal distribution has two parameters, you can't satisfactorily fit it to three constraints that don't naturally fit it. With extreme quantiles of 2.5 and 7.5, the mode is ~4, and there is not much you can do about it. Since the scale of the errors for a and b is much smaller than for c, one of them will be pretty much ignored during the optimization.

For a better fit, you can either pick a three-parameter distribution, for example the generalized gamma distribution (implemented in the VGAM package), or add a shift parameter to the lognormal (or gamma, ...) distribution.

As a last note, since the distribution you are looking for is clearly not symmetric, the average of the two given observations is not the right value for the mode. I would maximize the sum of the densities at 3.0 and 3.6 while maintaining the extreme quantiles at 2.5 and 7.5 - this is possible if you have three parameters.

-