# 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,ylim=c(0,4),p.beta,main="beta densities",type="l")

par(new=T)

beta1 <- 5
beta2 <- 1
p.beta <- dbeta(x,beta1,beta2)
plot(x,ylim=c(0,4),p.beta,main="beta densities",pch="+")

par(new=T)

beta1 <- 1
beta2 <- 5
p.beta <- dbeta(x,beta1,beta2)
plot(x,ylim=c(0,4),p.beta,main="beta densities",pch="*")

par(new=T)

beta1 <- 10
beta2 <- 10
p.beta <- dbeta(x,beta1,beta2)
plot(x,ylim=c(0,4),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))