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.
 A: 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.
A: 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.
