# beta(theta | a, b)
beta_pdf <- function(theta, a, b) {
  # Computes p(theta|a,b) = beta(theta|a,b)
  pdf = theta^(a - 1) * (1 - theta)^(b - 1)
  return(pdf)
}

theta = seq(0, 1, by = 0.01)
a = 2
b = 2
plot(theta,beta_pdf(theta,a,b))

beta_pde <- function(theta) {
  # Computes p(theta|a,b) = beta(theta|a,b)
  pdf <- seq(0, 100)
  for(i in 1:length(theta)){
  if(theta[i]==0.5){
    pdf[i] = 1
  }
  else{
    pdf[i]=0
  }
  }
  return(pdf)
}

plot(theta,beta_pde(theta))