I have a (probably) biased coin. I want to make a Bayesian inference on $\lambda=p(coin=\text{heads})$. I define my prior $p(\lambda)\sim Beta(\lambda; a=1,b=1)$. I toss the coin and update my posterior by simply increasing hyperparameters of the Beta distribution.

Now my question is if I know that $0.2\leq\lambda\leq0.7$, how can I make Bayesian inference? I tried to integrate Beta function for $[0.2, 0.7]$ to find normalizing constant. However I couldn't find out a compact form. I also don't know whether Beta distribution is an appropriate prior in this case.

  • $\begingroup$ What do you mean by "compact form"? $\endgroup$ – jbowman Oct 9 '18 at 15:24

You can perform the Bayesian updating of the parameters as you would with an unconstrained $\lambda$, then adjust the posterior to reflect the constraints by limiting the range of $\lambda$ to $[0.2,0.7]$ and renormalizing appropriately. In your case, as you realized, you'd have to be able to integrate the posterior Beta distribution to find the normalization constant. If, by good fortune or planning, you've picked a prior on $\lambda$ with integer parameters $a$ and $b$ (as you have done in the question), then the posterior will have integer parameters $a'$ and $b'$, and there is a closed form solution for the incomplete beta function:

$$I_x(a,b) = {1\over \text{B}(a,b)}\int_0^xp^a(1-p)^b\text{d}p=\sum_{j=a}^{a+b-1}{a+b-1\choose j}x^j(1-x)^{a+b-1-j}$$

... admittedly awkward if $a$ and/or $b$ are large.

A check of the formula, using R:

Ix <- function(x,a,b) {
  res <- 0
  m <- a + b - 1
  for (j in a:m) {
    res <- res + choose(m, j)*x^j*(1-x)^(m-j)

and comparison with the cumulative Beta distribution in the same language:

> Ix(0.4,3,5)
[1] 0.580096
> pbeta(0.4,3,5)
[1] 0.580096

In your case, the normalization of the posterior would involve dividing by $I_{0.7}(a',b')-I_{0.2}(a',b')$.

If, on the other hand, you don't have integer parameters, you will be forced to resort to using an infinite series expansion or numerical integration; the latter will probably work well enough given the smoothness of the functions.

Of course, languages such as R and Python / scikit have functions that will evaluate the CDF of the Beta distribution for you, so you can avoid the whole issue by simply using the available canned routines to find the normalization constant.

  • $\begingroup$ Well if I did it correctly, scipy gave some unreasonable results (NaN values where i thought it should not). My calculation was: $beta(a,b) - (betainc(0.2;a,b) * beta(a,b) + (beta(a,b)-betainc(0.7;a,b)*beta(a,b)))$ I think $betainc(x;a,b) * beta(a,b)$ is $I_x(a,b)$. $\endgroup$ – Alper Ahmetoglu Oct 9 '18 at 21:08
  • $\begingroup$ scipy.special.betainc actually gives you the normalized incomplete beta function, as in the formula above; this is equal to the cumulative density function of the Beta distribution, so all you need is betainc(a,b,0.7)-betainc(a,b,0.2) to get the normalization constant. Note the order of the parameters in the function call! $\endgroup$ – jbowman Oct 9 '18 at 22:09

' The standard way to estimate $\lambda$ would be to just use the estimator, would it not? $p(\mbox{heads}) = \frac{\mbox{# of heads}}{\mbox{# of flips}}$

"I toss the coin and update my posterior" Or are you trying to come up with an update to your estimate every single time you flip a coin, ie a recursive Bayesian estimation? You could do that by using the posterior at flip $n$ to be the prior distribution at flip $n+1$ and multiplying that by the likelihood to obtain the posterior at flip $n+1$. But the prior you give, a flat prior on $[0.2, 0.7]$ isn't conjugate to the Bernoulli likelihood, so I don't think you could make a nice updating chain..." `

Is there some particular reason why you want to do this recursively, updating your estimate with every new flip?

  • 2
    $\begingroup$ This doesn't answer the question, I'm sorry to say. $\endgroup$ – jbowman Oct 9 '18 at 16:52
  • $\begingroup$ Well, I do not need to update it everytime but I need posterior. I see that it is not conjugate to the Bernoulli likelihood but I thought that there should be some easy trick to do so, like integrating and finding some nice form. $\endgroup$ – Alper Ahmetoglu Oct 9 '18 at 21:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.