# a poisson regression using MCMC
# data from Dobson- An Introduction to Generalized Linear Models
y_scan()
2 3 6 7 8 9 10 12 15

x_scan()
-1 -1 0 0 0 0 1 1 1

n_length(y)
v0_100; v1_100;
sx_sum(x); sy_sum(y); sxy_sum(x*y)

pib0_function(b0){
  if(any(b0+b[2]*x < .1)) return(-9e9)
  h_-n*b0 - 1/(2*v0)*b0^2 + sum(y*log(b0+b[2]*x)) }
pib1_function(b1){
  if(any(b[1]+b1*x < .1)) return(-9e9)
  h_-sx*b1 - 1/(2*v1)*b1^2 + sum(y*log(b[1]+b1*x)) }

# initial values:
b_c(7,4)
# width of metropolis proposals:
r_c(1,1)

#number of iterations
N_500
burn_100
bs_matrix(nrow=N,ncol=length(b))

for(i in 1:N){
   # update b[1]
  bcan_runif(1,min=b[1]-r[1],max=b[1]+r[1])
  ldiff_pib0(bcan) - pib0(b[1])
  if(log(runif(1)) < ldiff) b[1]_bcan
   # update b[2]
  bcan_runif(1,min=b[2]-r[2],max=b[2]+r[2])
  ldiff_pib1(bcan) - pib1(b[2])
  if(log(runif(1)) < ldiff) b[2]_bcan
  bs[i,]_b
}

par(mfrow=c(2,2),mar=c(5,5,2,2),oma=c(1,1,3,1))
# show some realizations
plot(x,y); 
for(i in 1:10){
 rowval_sample(burn:N,1)
 abline(bs[rowval,1],bs[rowval,2])
}
mtext("posterior realizations",side=3,line=1)
# show the posterior mean fit
plot(x,y); abline(mean(bs[-(1:burn),1]),mean(bs[-(1:burn),2]))
mtext("posterior mean fit",side=3,line=1)

plot(bs[,1],bs[,2],type="p", xlab="intercept",ylab="slope") 
mtext("scatterplot of realizations",line=1)
matplot(1:N,bs,type="l"); mtext("time sequence of parameters",line=1)
mtext("Poisson regression with identity link",outer=T,side=3)


################# Now using the log link ########################
# a poisson regression using MCMC
# data from Dobson- An Introduction to Generalized Linear Models
y_scan()
0 1 2 3 1 4 9 18 23 31 20 25 37 45

x_log(1:14)


n_length(y)
v0_100; v1_100;
sx_sum(x); sy_sum(y); sxy_sum(x*y)

pib0_function(b0){
  h_sy*b0 - 1/(2*v0)*b0^2 - sum(exp(b0+b[2]*x)) }
pib1_function(b1){
  h_sxy*b1 - 1/(2*v1)*b1^2 - sum(exp(b[1]+b1*x)) }

# initial values:
b_c(-2,2)
# width of metropolis proposals:
r_c(.2,.1)

#number of iterations
N_5000
burn_100
bs2_matrix(nrow=N,ncol=length(b))

for(i in 1:N){
   # update b[1]
  bcan_runif(1,min=b[1]-r[1],max=b[1]+r[1])
  ldiff_pib0(bcan) - pib0(b[1])
  if(log(runif(1)) < ldiff) b[1]_bcan
   # update b[2]
  bcan_runif(1,min=b[2]-r[2],max=b[2]+r[2])
  ldiff_pib1(bcan) - pib1(b[2])
  if(log(runif(1)) < ldiff) b[2]_bcan
  bs2[i,]_b
}

par(mfrow=c(2,2),mar=c(5,5,2,2),oma=c(1,1,3,1))
# show some realizations
r_seq(0,3,length=50)
plot(x,y); 
for(i in 1:10){
 rowval_sample(burn:N,1)
 lines(x,exp(bs2[rowval,1]+bs2[rowval,2]*x))
}
mtext("posterior realizations",side=3,line=1)
# show the posterior mean fit
plot(x,y) 
lines(x,exp(mean(bs2[-(1:burn),1])+mean(bs2[-(1:burn),2])*x))
mtext("posterior mean fit",side=3,line=1)

plot(bs2[,1],bs2[,2],pch=".", xlab="intercept",ylab="slope") 
mtext("scatterplot of realizations",line=1)
matplot(1:N,bs2,type="l"); mtext("time sequence of parameters",line=1)
mtext("Poisson regression with log link",outer=T,side=3)