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.
I would like to take some slight shortcuts with regard to Romke's excellent solution while making use of the comfortable features of Mathematica's statistical framework.
$\theta|I \sim Beta(6,14)$
We can interpret this (conjugate) prior as having seen (or believing in) 18 Bernoulli-trials and having observed 5 successes and 13 failures before making the actual experiment (part of our background information $I$) and code this as:
priorDist = With[ { a = 6, b = 14 }, BetaDistribution[a, b] ];
prior = Function[ θ, PDF[ priorDist ] @ θ ];
Plot[ prior[θ], {θ, 0, 1}, Filling -> Axis, Axes-> {True,False} ]
$y|\theta, I \sim B(n, \theta)$
Given a series of $n = 1830$ Bernoulli-trials with $y = 420$ successes we will have a the Binomal PDF as a likelihood function:
likelihood = With[
{
n = 1830,
y = 420
},
Function[θ, Likelihood[ BinomialDistribution[n, θ], {y} ] ]
];
Plot[likelihood[θ], {θ, 0, 1}, Filling -> Axis, Axes -> {True, False}, PlotRange -> All]
As Romke pointed out, this function is a PDF depending on $\theta$ but not a probability distribution. So it does not integrate to 1:
NIntegrate[ likelihood[θ], {θ, 0, 1} ]
0.00054615
$p(\theta|y,I) \propto p(y|\theta,I) \cdot p(\theta|I)$
Bayes-Rule tells us, that the PDF of the posterior distribution equals the product of the likelihood function and the prior PDF up to a normalizing constant (often called the evidence), which Romke has used in his explicit answer:
$p(y|I) = \int_0^1 p(y|\theta,I) p(\theta|I)$
Being a bit more lazy we can make use of the fact that ProbabilityDistribution has an option Method -> "Normalize" which will take care to adjusting the PDF to meet the requirements of a probability distribution:
posteriorDistribution = ProbabilityDistribution[
prior[θ] likelihood[θ],
{θ, 0, 1},
Method -> "Normalize"
];
Plot[ PDF[ posteriorDistribution ][θ], {θ,0,1},
Filling->Axis, Axes->{True,False}, PlotRange->All]
After so many additional data, the posterior (of course) is very close to the likelihood. We can immediately get some summarizing statistics:
Through[{Mean, StandardDeviation}@posteriorDistribution] // N
{0.23027, 0.00978554}
RandomVariate will also work properly, so that we can sample from the posterior:
RandomVariate[ posteriorDistribution, 10 ]
{0.222989, 0.227517, 0.218535, 0.228404, 0.226061, 0.235919, 0.234529, 0.230528, 0.22929, 0.247735}
Since the prior and the posterior distributions here have the same functional form (conjugate prior), we may have immediately used the fact that:
$\theta|y,I \sim Beta(a + y, b + n - y)$
Equal[
PDF[ BetaDistribution[6 + 420, 14 + 1830 - 420], θ ],
PDF[posteriorDistribution, θ]
] // Reduce
True
You could also look at likelihood[θ] prior[θ]//Simplify to see that the posterior does have the form:
$constant \times \theta^\alpha (1 - \theta)^\beta$
Note that in general the normalizing constant (e.g. the evidence) can be obtained as an expectation:
$p(y|I) = \int p(y,\theta|I)d\theta = \int p(y|\theta,I) p(\theta|I)d\theta = \mathbb{E}[ p(y|\theta,I) ] = \mathbb{E}[\mathcal{L}(\theta)]$
posteriorPDF = With[
{ evidence = Expectation[ likelihood[θ], θ \[Distributed] priorDist ] },
Function[
θ,
prior[θ] likelihood[θ] / evidence
]
];
posteriorPDF[θ] == PDF[BetaDistribution[6 + 420, 14 + 1830 - 420], θ] // Reduce
True
Likelihood and Expectation. Doing so should be more general with regard to discrete or continuous distribution functions.
$\endgroup$
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])
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}
I created a snippet that generates the posterior predictive distribution via the above method:
BayesRule[priorDist_?DistributionParameterQ,
sampleDist_?DistributionParameterQ, params_List, data_List] :=
Module[
{bounds, boundReg, likelihood, prior, posterior},
boundReg =
RegionBounds[
ImplicitRegion[DistributionParameterAssumptions[sampleDist],
params], "Sufficient"];
bounds = Flatten /@ Transpose[{params, boundReg}];
prior = PDF[priorDist];
likelihood =
Function[Evaluate[params], Likelihood[sampleDist, data]];
posterior =
ProbabilityDistribution[(prior @@ params)*(likelihood @@ params),
Sequence @@ bounds, Method -> "Normalize"];
Return[ParameterMixtureDistribution[sampleDist,
If[Length[params] == 1, params[[1]], params] \[Distributed]
posterior]]
]
For example, to solve the problem in the original question:
BayesRule[BetaDistribution[6, 14],
BinomialDistribution[1830, \[Theta]], {\[Theta]}, {420}]
Which creates a distribution that is equivalent to the beta-binomial distribution
