help(dbeta)to get help about the pdf for the Beta distribution.
Notice that the help message gives help on several functions related to the Beta. This bonus information is how R creates help messages for any of its built in probability distributions. Click Here to see a list of available distributions.
beta.plot <- function(r,s) {
x <- seq(0,1,length=501)
y <- dbeta(x,r,s)
plot(x,y,type='l',main=paste("r =", r, ",", "s =", s))
}
Enter that function into R. Then run the following code to see a collection of beta pdf's.
m <- 4
n <- 4
par(mfrow=c(m,n))
r.vect <- c(.75, 2, 5, 100)
s.vect <- r.vect
for (i in 1:m) {
for (j in 1:n) {
r <- r.vect[i]
s <- s.vect[j]
beta.plot(r,s)
}
}
par(mfrow=c(1,1))
Thumbtacks : Empirically, I dropped a tack 50 times and got 35 "point-ups". So n=50 and k=35. The code below plots the appropriate posterior beta distribution using different beta priors.
k <- 35 n <- 50 r <- 1; s <- r # uniform prior par(mfrow=c(2,1)) beta.plot(r,s) beta.plot(k+r,n-k+s) qbeta(.025,k+r,n-k+s) # lower limit of posterior interval qbeta(.975,k+r,n-k+s) # upper limit of posterior interval (k+r)/(n+r+s) # posterior mean, aka "squared-error loss Bayes estimate r <- 10; s <- r # symmetric about .5; sd=.11 par(mfrow=c(2,1)) beta.plot(r,s) beta.plot(k+r,n-k+s) qbeta(.025,k+r,n-k+s) # lower limit of posterior interval qbeta(.975,k+r,n-k+s) # upper limit of posterior interval (k+r)/(n+r+s) # posterior mean, aka "squared-error loss Bayes estimate r <- 100; s <- r # symmetric about .5, sd=.035 par(mfrow=c(2,1)) beta.plot(r,s) beta.plot(k+r,n-k+s) qbeta(.025,k+r,n-k+s) # lower limit of posterior interval qbeta(.975,k+r,n-k+s) # upper limit of posterior interval (k+r)/(n+r+s) # posterior mean, aka "squared-error loss Bayes estimate r <- 1000; s <- r # symmetric about .5, sd=.011 par(mfrow=c(2,1)) beta.plot(r,s) beta.plot(k+r,n-k+s) qbeta(.025,k+r,n-k+s) # lower limit of posterior interval qbeta(.975,k+r,n-k+s) # upper limit of posterior interval (k+r)/(n+r+s) # posterior mean, aka "squared-error loss Bayes estimate Here is a summary table of the results:Beta Prior Estimate Lower Limit Upper Limit r=s=1 .692 .562 .809 r=s=10 .643 .528 .750 r=s=100 .540 .478 .601 r=s=1000 .505 .483 .527 Wald int .700 .573 .827 Plus-4 .685 .561 .809