# When using the beta distribution as a prior distribution for binomial, why won't the distribution results match with the calculated probability?

Let's say I have 1 success in 4 bernoulli trials, and I wish to plot the distribution of the parameter $p$ of the corresponding binomial distribution. I'm using R.

The probability of seeing 1 sucess and 3 failures in 4 tests for $p=0.25$ is, for these parameters:

> n <- 4
> p <- 0.25
> dbinom(1, n, p)
[1] 0.421875


To get the distribution for the parameter, I use a beta distribution $Beta(k+1, n-k+1)$. But when I try to calculate the value for $p=0.25$, I get a different result:

> k <- 1
> dbeta(p, k+1, n-k+1)
> [1] 2.109375


I tried to divide that by the function beta $B(k+1, n-k+1)$, but that did not work either:

> dbeta(p, k+1, n-k+1)/beta(k+1, n-k+1)
> [1] 42.1875


That's exactly 100 times the result I expected. I am certainly missing something - maybe something basic, I'm no expert.

Another thing that intrigues me: dbeta region sums to 1 without dividing by the beta function.

You are misunderstanding the posterior distribution. The height of the posterior is mainly useful for plotting, the value at a single plot is not interpretable. The area under the curve for a given range represents the posterior probability that the parameter is in that range, so you need to integrate over a small range to get the probability.

The probability you computed using dbinom is the probability of getting 1 success given that the parameter is 0.25, which is very different from the posterior on the parameter given the data. Since they are completely different things, there is no reason to expect them to match.

• Thanks! I understand that the density function needs to be integrated over the range I want. But I am curious: what do I have if I plot the probability of getting one success in 4 trials for a lot of $p$ values? A PMF? – Penz Aug 14 '11 at 20:47
• That is the liklihood function. If you scale it to have total area of 1 then it could be the posterior from using a Beta(0,0) (improper) prior. – Greg Snow Aug 14 '11 at 23:25

You will get a more intuitive result by looking at the posterior mode, which is given by $p=\frac{1}{4}$.

You say at the start that you want to plot the distribution, but then you start calculating things, rather than plotting them. It would be more straight forward to create the plot:

k <- 1
n <- 4
p <- 0:100/100
density <- dbeta(p, k+1,n-k+1)
plot(denisty~p)


Note that the proper, mathematically precise way to write the density is:

$$Pr\big(p\in(p_{0},p_{0}+dp)|nkI\big)=\frac{p_{0}^{k}(1-p_{0})^{n-k}}{B(k+1,n-k+1)}dp$$

The R function "dbeta" gives you the the fraction, but you need to multiply by the differential $dp$ in order to get a probability. The distribution function or cummulative density function is what is a probability (and will be between $0$ and $1$), and this is given by the "pbeta" R function. So you could plot this as well (using p, n, and k from above code)

CDF <- pbeta(p,k+1,n-k+1)
plot(CDF~p)


The CDF is the probability that $0<p<p_0$

dbeta already includes the beta function in the denominator, and you can check this manually for we have $p_{0}=0.25=\frac{1}{4}$ so $p_{0}^{1}(1-p_{0})^{3}=\frac{3^{3}}{4^{4}}=\frac{27}{256}$. The beta function is given by $B(2,4)=\frac{1!3!}{5!}=\frac{1}{20}$. So we get the density as $\frac{27\times 20}{256}=\frac{135}{64}=2.109375$

On your last point, the dbeta region should sum to $1$, because we have a normalised probability density. This is basically $\int_{0}^{1}f(p)dp=1$, and we also have $f(p)\geq 0$, so the integral is the area between the density, plotted on the y-axis and the x-axis.