Bayesian Inference with Continuous prior distribution

Question

Assume that we have a prior distribution for the probability of success as follows.

p = BetaDistribution[6, 14]

which we can graph as this

Plot[Evaluate[PDF[p, k]], {k, 0, 1}]

I would like to use Mathematica to calculate the posterior distribution given that we take a sample and we get 420 successes out of a sample size of 1830.

what would be to calculate it and graph it?

Regards.

Answer 1

The posterior distribution for p is a Beta [426, 1424]. You can graph it as you did yourself.

Here is how you can compute in Mathematica.

In[1]:= bin = Binomial[n, x]*p^x*(1 - p)^(n - x)

Out[1]= (1 - p)^(n - x) p^x Binomial[n, x]

In[2]:= prior = p^(a - 1)*(1 - p)^(b - 1)/Beta[a, b]

Out[2]= ((1 - p)^(-1 + b) p^(-1 + a))/Beta[a, b]

In[3]:= den = 
 Integrate[bin*prior, {p, 0, 1}, 
  Assumptions -> {a > 0, b > 0, x >= 0, n >= x}]

Out[3]= (Binomial[n, x] Gamma[b + n - x] Gamma[a + x])/(
Beta[a, b] Gamma[a + b + n])

In[4]:= bin*prior/den // FullSimplify

Out[4]= ((1 - p)^(-1 + b + n - x) p^(-1 + a + x) Gamma[a + b + n])/(
Gamma[b + n - x] Gamma[a + x])

Answer 2

I would like to extend asim's solution, which is a function, to a probability distribution.

ClearAll[evidence, likelihood, prior, posterior];
likelihood = p^x*(1 - p)^(n - x) (* is always unnormalized ! *)
prior = p^(a - 1)*(1 - p)^(b - 1)/Beta[a, b] (* always normalized *)
evidence = Integrate[likelihood*prior, {p, 0, 1},
Assumptions -> {a > 0, b > 0, x >= 0, n >= x}] (* normalisation of posterior *)
posterior = prior*likelihood/evidence // FullSimplify

(* Out[]= *)
(1 - p)^(n - x) p^x
((1 - p)^(-1 + b) p^(-1 + a))/Beta[a, b]
(Gamma[b + n - x] Gamma[a + x])/(Beta[a, b] Gamma[a + b + n])
((1 - p)^(-1 + b + n - x) p^(-1 + a + x) Gamma[a + b + n])/(Gamma[b + n - x] Gamma[a + x])

Nothing new so far. Now for the probability distribution:

ClearAll[posteriorDistribution];
posteriorDistribution[a_, b_, x_, n_] = ProbabilityDistribution[posterior, {p, 0, 1}, 
Assumptions -> {a > 0, b > 0, x >= 0, n >= x}];
PDF[posteriorDistribution[6, 14, 420, 1830], 0.21] (* test at p=0.21 *)

4.6369

Check this answer by asim's solution in BetaDistribution[] form:

PDF[posteriorDistribution[6, 14, 420, 1830], 0.21] == PDF[BetaDistribution[426, 1424], 0.21]

True

This probability distribution can be plotted by:

Plot[PDF[posteriorDistribution[6, 14, 420, 1830], p], {p, 0, 1}, PlotRange -> All]

Several statistical properties work fine:

Mean[posteriorDistribution[6, 14, 420, 1830]] // N
StandardDeviation[posteriorDistribution[6, 14, 420, 1830]] // N

0.23027
0.00978554

But RandomVariate[] does not:

RandomVariate[posteriorDistribution[6, 14, 420, 1830], 5]

RandomVariate[
 ProbabilityDistribution[
 25100175092366680639705706981674366907407305357564535864428789165880\
7774460815201505720742948687064397263210935672027225834571883138571552\
9435073323698076733502611129508985779654636172308803925192878804616791\
4163296382700546574839309987790038929990089862577986787654011552037119\
0043901814728663569787527433974843693348538658624431557074602596292870\
1098499374967102449446582880605618133093761170547891966862458477624898\
28252384878666048 (1 - \[FormalX])^1423 \[FormalX]^425, {\[FormalX], 
0, 1}, Assumptions -> {True, True, True, True}], 5]

Anyone here who knows what is going wrong?

Of course, the BetaDistribution[] works fine:

RandomVariate[BetaDistribution[426, 1424], 5]

{0.224794, 0.246947, 0.225743, 0.241723, 0.219218}