## October 10, Bayesian Estimation 
rm(list = ls())

###################################################################################
## Poisson Example
salmon.count <- rpois(31, 10)  ## generate fake data for a month
y <- salmon.count

## we know the mle is the mean
mle.lambda <- mean(y)
mle.lambda

## set our prior beliefs about the number of fish per day

## let's say we are not knowledgable about fish and think the distribution can vary quite a bit
par(mfrow=c(2,1))
alpha.diff <- 3
beta.diff <- 4
hist(rgamma(1000, shape=alpha.diff, scale=beta.diff), col="yellow", main="Prior Beliefs about the Number of Fish")

## let's draw 1000 samples from the posterior distribution
alpha.star <- sum(y)+alpha.diff
beta.star <- beta.diff/(beta.diff*length(y)+1)

post.lambda <- rgamma(1000, shape=alpha.star, scale=beta.star)
hist(post.lambda, col="yellow", main="Posterior Distribution")
abline(v=mean(post.lambda), col="red", lwd=3)
abline(v=mean(mle.lambda), col="blue", lwd=3)

## Now we are fish experts!!  We know very well the number of fish over the ladder every day
## mean = 5, var = 1
## thus
alpha.know <- 25
beta.know <- 1/5
par(mfrow=c(2,1))
hist(rgamma(1000, shape=alpha.know, scale=beta.know), col="yellow", main="Prior Beliefs about the Number of Fish")


## let's draw 1000 samples from the posterior distribution
alpha.star <- sum(y)+alpha.know
beta.star <- beta.know/(beta.know*length(y)+1)

post.lambda <- rgamma(1000, shape=alpha.star, scale=beta.star)
hist(post.lambda, col="yellow", main="Posterior Distribution")
abline(v=mean(post.lambda), col="red", lwd=3)
abline(v=mean(mle.lambda), col="blue", lwd=3)
mean.post.lambda.31 <- mean(post.lambda)


## Now suppose we have a lot more data, maybe a year's worth
## Poisson Example
par(mfrow=c(1,1))
salmon.count <- rpois(365, 10)  ## generate fake data for a month
y <- salmon.count

## let's draw 1000 samples from the posterior distribution
alpha.star <- sum(y)+alpha.know
beta.star <- beta.know/(beta.know*length(y)+1)

post.lambda <- rgamma(1000, shape=alpha.star, scale=beta.star)
hist(post.lambda, col="yellow", main="Posterior Distribution", xlim=c(8.5,11))
abline(v=mean(post.lambda), col="red", lwd=3)
abline(v=mean(mle.lambda), col="blue", lwd=3)
abline(v=mean.post.lambda.31, col="green", lwd=3)