# Help in Bayesian Bernoulli-Beta Model (solution verification)

I have some doubts that -for many of you- might seem basic and I am basically looking for some guidance verifying my attempt of solution.

The Problem

Let $$p(x\mid\theta)=Bernoulli(\theta)$$, and let's suppose the Prior Distribution is given by $$p(\theta)=Beta(\theta \mid \alpha,\beta)$$. ($$\theta$$ is unknown, and $$\alpha, \beta$$ are known).

Assume independet observations and obtain: 1) The Posterior Distribution of $$\theta$$, 2)The Prior Predictive Distribution, and 3)The Posterior Predictive Distribution.

Attempt of Solution

Let $$I=\{1,2,...,n\}$$, $$\mathrm{X}=\{X_i\}_{i \in I}$$ the independent observatrions above mentioned, and $$\mathbb{I}$$ be the Indicator function.

1)

Given the information of the excercise,

$$p(\theta \mid \mathrm{x})=\frac{[\prod_{i \in I}\theta^{x_i}(1-\theta^{1-{x_i}})\mathbb{I}_{\{0,1\}}(x_i)][\theta^{\alpha-1}(1-\theta)^{1-\beta} \mathbb{I}_{[0,1]}(\theta)]}{\int_\Theta [\prod_{i \in I}\theta^{x_i}(1-\theta^{1-{x_i}})\mathbb{I}_{\{0,1\}}(x_i)][\theta^{\alpha-1}(1-\theta)^{1-\beta} \mathbb{I}_{[0,1]}(\theta)]d\theta}$$

$$p(\theta \mid \mathrm{x})=\frac{\theta^{\alpha + r-1}(1-\theta)^{\beta+n-r-1}\mathbb{I}_{\{0,1\}}(x_i)\mathbb{I}_{[0,1]}(\theta)}{\int_0^1\theta^{\alpha + r-1}(1-\theta)^{\beta+n-r-1}\mathbb{I}_{\{0,1\}}(x_i)d\theta}$$

$$p(\theta \mid \mathrm{x})=Beta(\theta \mid \alpha+r,\beta+n-r),$$ $$r:=\sum_{i \in I}x_i$$.

2) The Prior Predictive Distribution $$p(x)$$ is given then by

$$p(x)=\frac{1}{Beta(\theta \mid \alpha,\beta)}\int_0^1\theta^{\alpha+x-1}(1-\theta)^{\beta-x}\mathbb{I}_{\{0,1\}}(x)d\theta$$

$$p(x)=\frac{1}{Beta(\theta \mid \alpha,\beta)}\int_0^1(\theta^{\alpha+x-1}-\theta^{\alpha+\beta-1})\mathbb{I}_{\{0,1\}}(x)d\theta$$

$$p(x)=\frac{1}{Beta(\theta \mid \alpha,\beta)}(\frac{1}{\alpha+x}-\frac{1}{\alpha+\beta})\mathbb{I}_{\{0,1\}}(x)$$

3) The Posterior Predictive Distribution $$p(x \mid \mathrm{x})$$ is given then by

$$p(x \mid \mathrm{x})=\frac{1}{Beta(\theta \mid \alpha+r,\beta+n-r)}\int_0^1\theta^{\alpha+x+r-1}(1-\theta)^{\beta+n-x-r}\mathbb{I}_{\{0,1\}}(x)d\theta$$

$$p(x \mid \mathrm{x})=\frac{1}{Beta(\theta \mid \alpha+r,\beta+n-r)}\int_0^1(\theta^{\alpha+x+r-1}-\theta^{\alpha+\beta+n-1})\mathbb{I}_{\{0,1\}}(x)d\theta$$

$$p(x \mid \mathrm{x})=\frac{1}{Beta(\theta \mid \alpha+r,\beta+n-r)}(\frac{1}{\alpha+x+r}-\frac{1}{\alpha+\beta+n})\mathbb{I}_{\{0,1\}}(x)$$

I have not verified every line of your Question, but it seems your initial formulation for each part is on the right track. In this particular situation, you may be able to save some work by noticing that the beta prior and the binomial likelihood are 'conjugate', as discussed below. [I will not comment on the posterior predictive distribution because you have not defined it or explained the meaning of $$r.]$$

Prior. As a specific example, suppose you are trying to interpret a public opinion poll to judge which of two candidates for mayor of a large city is currently ahead. Many American elections result in something near a 50:50 vote count. So it might be reasonable to use the prior distribution $$\mathsf{Beta}(5, 5),$$ which puts nearly half of its probability between 40% and 60% of the vote. [Computation in R.]

diff(pbeta(c(.4,.6), 5, 5))
[1] 0.4668646


If $$\theta$$ is the percentage in favor of Candidate A, then we might write this as $$p(\theta) = K\theta^{5-1}(1-\theta)^{5-1},$$ where $$K$$ is the norming constant (constant of integration). Sometimes we write this as $$p(\theta) \propto \theta^{5-1}(1-\theta)^{5-1},$$ where the symbol $$\propto$$ (read 'proportional to') acknowledges than the factor $$K$$ is absent.

Likelihood. Then suppose results from a reliable pollster based on interviewing $$n = 700$$, reports $$x = 372$$ in favor of Candidate A. Then the binomial likelihood function for the poll would be $$p(x|\theta)\propto \theta^{372}(1-\theta)^{700-372}.$$

Posterior. Then by Bayes' Theorem we have that the posterior distribution $$p(\theta|x) = p(\theta) \times p(x|\theta).$$ In our case this gives $$p(\theta|x) = \theta^{5-1}(1-\theta)^{5-1} \times \theta^{372}(1-\theta)^{700-372} \propto \theta^{377-1}(1-\theta)^{333-1},$$ where we recognize the right-hand side as proportional to the density function of $$\mathsf{Beta}(377, 333).$$

Conjugacy. Because the beta and binomial distributions are 'conjugate' (mathematically compatible), we do not need to do further computation to find the norming constant for the posterior distribution.

Posterior probability interval. A 95% posterior probability interval for the population proportion in favor of Candidate A would be $$(0.494, 0.568).$$ Disappointingly, this interval includes 0.5, so the candidate should maintain his/her campaign, but chances of success seem good.

qbeta(c(.025, .975), 377, 333)
[1] 0.4942327 0.5675738


Note: For the poll discussed above, a frequentist Agresti-Coull 95% confidence interval is $$(0.4944, 0.5681),$$ numerically about the same as the Bayesian posterior (or 'credible') interval, but subject to somewhat different interpretations.

The length of the standard (probability-symmetric) interval above, which puts probability $$0.025$$ in each tail, is:

std.ci = qbeta(c(.025,.975), a1, a2); diff(std.ci)
[1] 0.07334116


Your posterior distribution with shape parameters 377 and 333, which have a ratio of about 1, is very nearly symmetrical, so the shortest interval will not be much different from the standard interval:

a1 = 377; a2 = 333; cred = .95
ltp = seq(.0001, 1-cred, by=.00001) # vector of lower tail probs
ll = qbeta(ltp, a1,a2); ul=qbeta(cred+ltp, a1,a2)
len = ul-ll
ltp.short = ltp[len==min(len)];  ltp.short
[1] 0.02518    # not much different from 0.025
short.ci = qbeta(c(ltp.short, cred+ltp.short),a1,a2)
diff(short.ci)
[1] 0.07334081 # not much shorter than standard int


As @Xi'an comments, this shortest interval is called the HPD interval, for 'highest probability density', for reasons covered in various comments.

• Thanks a lot for such a clear explanation and. This is more than I expected. Thank you. – Israel Barquín Oct 25 '19 at 21:07
• I have one more question: is there a general non-numerical method to find a minimal-length interval in the Bayesian approach? I have one problem in which I have to find such region given a Beta Posterior Distribution for my unknown parameter and I am having such a hard time with that, stats.stackexchange.com/users/85665/bruceet. – Israel Barquín Oct 27 '19 at 2:05
• There are various methods for finding the shortest posterior interval. I usually write a small program that cuts various probabilities from the left tail, incrementing by 0.1% or 0.01%, and do a search for the shortest. (Maybe something more elegant is programmed into one of the R libraries.) – BruceET Oct 27 '19 at 2:12
• That is what I was worried about. In a couple of days I will take an exam about all of this. I was given a guide in which some problems require me to find such intervals. I won't have computational instruments during this test. – Israel Barquín Oct 27 '19 at 2:19
• Seem to recall that the height of beta posterior is the same at the quantiles that give the endpoints of the shortest posterior interval. – BruceET Oct 27 '19 at 2:21

The answer to (2) is incorrect since $$\theta^{\alpha+x-1}(1-\theta)^{\beta-x}\ne\theta^{\alpha+x-1}-\theta^{\alpha+\beta-1}$$ Similarly, the answer to (3) is incorrect due to the same mistake.

• Oh, yes. What an embarrassing mistake. I will fix it, Thank you. – Israel Barquín Oct 27 '19 at 1:53
• Thanks for mention of errors (+1). // Still waiting for OP to fix this. – BruceET Oct 27 '19 at 16:36