# com26a of september 2021
# with Nils Exercise 1, stk 4021

# gamma prior, then poisson data, and updating, point by point 
# basic computation: 
# theta is gamma(a,b) implies 
# theta | y_1, ..., y_n is gamma(a + z_n, b + n), with z_n = sum of data 

# gamma prior (a,b): 
aa <- 0.1
bb <- 0.1 

# can try an overconfident prior, too sure of itself: 
# e.g. mean = 10, small variance: 
# aa <- 100*10.0
# bb <- 100 

yy <- c(6, 8, 7, 6, 7, 4, 11, 8, 6, 3)
# i did this from
# yy <- rpois(25,6.666)
nn <- length(yy) 

thetaval <- seq(1.5,12.0,by=0.01)

# the so-called silly prior, clashing with data, 
# which have to learn, step by step: 
# aa <- cc*theta0 
# bb <- cc 

# aa <- 0.50*50
# bb <- 0.50

# computing lots of p(\theta | y_1,...,y_j) 
curves <- 0*thetaval %*% t(1:nn)

for (j in 1:nn)
{ 
curves[ ,j] <- dgamma(thetaval,aa+cumsum(yy)[j],bb+j)
}

colours <- sample(2:4,nn,replace=T)
linetypes <- sample(2:4,nn,replace=T) 
colours[1] <- 1 
colours[nn] <- 1 
linetypes[1] <- 1 
linetypes[nn] <- 1 

par(bty="l")
par(cex.lab=1.4142)

matplot(thetaval,curves,type="l",lwd=2.5,
   lty=linetypes,col=colours, 
   xlab="theta",ylab="posteriors, sequentially")

## then examining results from a "stupid prior", 
#  with mean 40 and stdev 10,
#  whereas thetatrue = 10: 

theta0 <- 40
cc <- theta0/100

# gamma(cc*theta0, cc): variance is 100, stdev is 10 
# gamma(16, 0.4) has mean 40 and variance 16/0.16 = 100 

thetatrue <- 10 
yy <- rpois(100,thetatrue) 

thetaval <- seq(5,60,by=0.1)

# after j data points: 
j <- 10
aprime <- cc*theta0 + j*mean(yy[1:j])
bprime <- cc + j 

priorval <- dgamma(thetaval,cc*theta0,cc)
posteriorval <- dgamma(thetaval,aprime,bprime)

matplot(thetaval,cbind(posteriorval,priorval),
   type="l",lwd=2.0,col=c(1,2))