# Plot the density function of the beta distribution for different values of
# beta1 and beta2.

x = (1:99)/100

beta1 = 1
beta2 = 1
p.beta = dbeta(x,beta1,beta2)
plot(x,p.beta,main="beta densities",type="l")

par(new=T)

beta1 = 5
beta2 = 1
p.beta = dbeta(x,beta1,beta2)
plot(x,p.beta,main="beta densities",pch="+")

par(new=T)

beta1 = 1
beta2 = 5
p.beta = dbeta(x,beta1,beta2)
plot(x,p.beta,main="beta densities",pch="*")

par(new=T)

beta1 = 10
beta2 = 10
p.beta = dbeta(x,beta1,beta2)
plot(x,p.beta,main="beta densities",pch="-")

# Take a random sample from a beta.

beta1 = 5
beta2 = 1
n = 1000
y = rbeta(n,beta1,beta2)
hist(y,main="histogram of a random sample from a beta distribution",xlim=c(0,1))