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

Question

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!

Answer 1

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?

Answer 2

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.

Answer 3

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.

Answer 4

You could also try the triangular distribution. To fit this, you basically specify a lower bound (this would be X=2), an upper bound (this would be X=8), and a "most likely" value. The wikepedia page http://en.wikipedia.org/wiki/Triangular_distribution has more info on this distribution. If there is not much faith in the "most likely" value (as it appears to be, prior to observing any data), it may be a good idea to place a non-informative prior distribution on it, and then use the two data points to estimate this value. One good one is the jeffrey's prior, which for this problem would be p(c) = 1/(pi*sqrt((c-2)*(c-8))), where "c" is the "most likely value" (consistent with the wikipedia notation).

Given this prior, you can work out the posterior distribution of c analytically, or via simulation. The analytic form of the likelihood is not particularly nice, so simulation seems to be more attractive. This example is particularly well suited to rejection sampling (see wiki page for a general description of rejection sampling), because the maximised likelihood is 1/3^n regardless of the value of c, which provides the "upper bound". So you generate a "candidate" from the jeffrey's prior (call it c_i), and then evaluate the likelihood at this candidate L(x1,..,xn|c_i), and divide by the maximised likelihood, to give (3^n)*L(x1,..,xn|c_i). You then generate a U(0,1) random variable, and if u is less than (3^n)*L(x1,..,xn|c_i), then accept c_i as a posterior sampled value, otherwise throw away c_i and start again. Repeat this process until you have enough accepted samples (100, 500, 1,000, or more depending on how accurate you want). Then just take the sample average of whatever function of c you are interested in (the likelihood of a new observation is an obvious candidate for your application).

An alternative to accept-reject is to use the value of the likelihood as a weight (and don't generate the u), and then proceed with taking weighted averages using all candidates, rather than un-weighted averages with the accepted candidates