2

I am creating a model that utilize the output probability distribution of another model, as an input. The model outputs a probability of an event occurring per year, the distribution is a probability between (between 0 and 1) and right skewed, it is defined by Mean = 0.019, Mode = 0.009, 90% confidence interval: 5th percentile = 0.0002, 95th percentile = 0.07. I want to use a beta distribution so I can perform a Bayesian update using a binomial likelihood for the update.

However, I am uncertain of 1) if it is appropriate to use a beta distribution (as opposed to a lognormal) and how to prove if it is appropriate or not, and 2) how to determine the parameters of the beta distribution from the available information.

In relation to 2) I have looked at several articles on cross validate which seem related but I cannot quite work out how to implement: Calculate the confidence interval for the mean of a beta distribution Calculating the parameters of a Beta distribution using the mean and variance

This one seems the closest and easiest to implement: Selecting alpha and beta parameters for a beta distribution, based on a mode and a 95% credible interval

In the model I am creating it is possible to define a reparametrized beta distribution according to Mean and standard deviation, so I thought that I would be able to solve equations with the above mentioned values to obtain the standard deviation and define the beta distribution but now I am not sure if this is possible.

noideaboutlogs
  • 31
  • 3
  • 2
    How "appropriate" an approximation is depends on how you assess the error in the approximation. For instance, if you use a log-Normal approximation for a distribution over $(0,1)$, how would you penalise the fact that this distribution produces realisations over $(0,\infty)$? To answer the question one would need this measure of error, for instance a Kellback-Leibler or a Wasserstein distance. – Xi'an Feb 03 '22 at 08:24
  • 1
    What exactly do you want to do? – Tim Feb 03 '22 at 08:40

2 Answers2

2

A beta distribution has support on $[0,1]$ while a log-normal distribution has support on $(0,\infty)$ so it is likely that at least one would be inappropriate for your problem, possibly both. In particular the log-normal might look a strange choice for modelling a probability.

You have four pieces of information, which is in a sense too much information for finding two parameters of a beta distribution. So you could only find an approximation, and optimising it would depend on how important the differences are.

For example $\alpha=1.5,\beta = 75$ would give

  • a mean of about $0.0196$, close to your figure
  • a mode of about $0.00671$, lower than your figure
  • a 5th percentile of about $0.00234$, more than ten times your figure
  • a 95th percentile of about $0.0506$, lower than your figure

and you would need to judge whether these (or those for some other choices of $\alpha$ and $\beta$) were good enough for your purposes. If not, then using a beta distribution may not be what you want

Henry
  • 30,848
  • 1
  • 63
  • 107
2

A beta distribution or (if a single beta does not do an adequate job) a mixture of beta distributions is quite a sensible choice for quantifying beliefs/knowledge/information about the likely values of a probability (while a log-normal is a less natural choice, because it would consider it possible for the probability to be above 1).

As discussed e.g. by Diaconis and Ylvisaker (1985) a mixture of conjugate priors can approximate any (practically relevant) distribution you want arbitrarily well. The nice thing is that you can even still do conjugate updating with such a mixture: You just do conjugate updating of each of the mixture components and the only tricky bit is how to update the weights of each mixture component. That gets done by multiplying the old weight with the probability mass function of the posterior predictive distribution (for a binomial outcome given the posterior mixture component for the parameter after conjugate updating), which gets you a beta-binomial pmf. Then you normalize the weights. Et voilà, your posterior.

E.g. the RBesT R package provides convenient facilities for doing this:

library(RBesT)
library(ggplot2)
library(patchwork)

# Sample some data from a normal and then use the 
# inverse-logit function to get probabilities.
samples = plogis(rnorm(n=10000, mean=-1, sd=0.5))
hist(samples, breaks=50)

# Find a 3-component mixture to approximate the samples
emfit1 = mixfit(samples, type = "beta", Nc=3)
# For plotting purposes we turn this into a mixture object
# if we wanted to continue the analysis we could just use the emfit object
mixture1 = mixbeta(comp1=emfit1[1:3,1], 
                   comp2=emfit1[1:3,2], 
                   comp3=emfit1[1:3,3], param="ab")
p1 = plot(mixture1, fun=dmix, prob=1)

# Posterior updating with that mixture prior
# and observed data of 10 successes out of 100 tries.
mixture2 = postmix(mixture1, n=100, r=10)
print(mixture2)
p2 = plot(mixture2, fun=dmix, prob=1-1e-13)

(p1 + ggtitle("Prior mixture)) + (p2 + ggtitle("Posterior mixture"))

enter image description here

As you can see, in the example the mixture component that fit the data best got the most weight.

References

Diaconis, P. and Ylvisaker, D. (1985). Quantifying prior opinion. In Bernardo, J. M., DeGroot, M. H., Lindley, D. V., and Smith, A. F. M., editors, Bayesian statistics 2: Proceedings of the second Valencia International Meeting, September 6-10, 1983, pages 133–156, North-Holland, Amsterdam.)

Björn
  • 21,227
  • 2
  • 26
  • 65