library(LearnBayes)

# Q.1

i <- 1:18
y <- c(15,11,14,17,5,11,10,4,8,10,7,9,11,3,6,1,1,4)
pois.data <- data.frame(i  = i,y = y)


# Aleph

plot(i,log(y))
abline(lm(log(y) ~ i))
summary(lm(log(y) ~ i))

pois.loglin <- function (theta, data) 
{
    i.vec <- data[, 1]
    n.vec <- data[, 2]

    log.lambda <- theta[1] + theta[2] * i.vec
    d.pois <- sum( n.vec * log.lambda  - exp(log.lambda))	
    return(d.pois)
}

points(2.91,-0.111,col=2,pch="+",cex=1.5)


# Bet 

mycontour(pois.loglin,c(1.5,4,-.2,.1),pois.data)
title(xlab="Intercept", ylab="slope")


(fit <- laplace(pois.loglin,c(3,-.1),pois.data))

npar <- list(m=fit$mode,v=fit$var)
mycontour(lbinorm,c(1.5,4,-.2,.1),npar,add=T,col=3)
points(fit$mode[1],fit$mode[2],col=3,pch="+",cex=1.5)

fit$mode[2] ; 1.96 * sqrt(fit$var[2,2]) ; fit$mode[2] - 1.96 * sqrt(fit$var[2,2])


#  Haashara: what eactly is the laplace approximation plotted in green ?

library(mvtnorm)

pois.loglin(fit$mode,pois.data)
dmvnorm(fit$mode,fit$mode,fit$var,log=T)

pois.loglin(fit$mode,pois.data) - dmvnorm(fit$mode,fit$mode,fit$var,log=T)
fit$int


#  Gimmel 


# 1. find poposal distribution

tpar    <- list(m=fit$mode,var=fit$var,df=4)
datapar <- list(data=pois.data,par=tpar)

pois.t.comp <- function(theta,data)
{
	data <- datapar$data
	tpar <- datapar$par
	d    <- pois.loglin(theta,data)-dmt(theta,mean=c(tpar$m),
							  S=tpar$var,df=tpar$df,log=TRUE)
return(d)
}

mycontour.DY <- 
function (logf, limits, data) 
{
    ng = 50
    x0 = seq(limits[1], limits[2], len = ng)
    y0 = seq(limits[3], limits[4], len = ng)
    X = outer(x0, rep(1, ng))
    Y = outer(rep(1, ng), y0)
    n2 = ng^2
    Z = logf(cbind(X[1:n2], Y[1:n2]), data)
    Z = Z - max(Z)
    Z = matrix(Z, c(ng, ng))
    contour(x0, y0, Z,levels = seq(-30,30,by = 3), lwd = 2)
}

mycontour.DY(pois.t.comp,c(-0,6,-2,2),pois.data)


# 2. Sample initial theta values from proposal distribution

theta.t.sample <- rmt(10^5, mean = c(tpar$m), S = tpar$var, df = tpar$df)
plot(theta.t.sample[,1],theta.t.sample[,2])
add.contour(pois.loglin,c(1.5,4,-.2,.1),pois.data)

d.vec <- rep(NA,10^5)
for(i in 1:10^5)
	d.vec[i] <- pois.loglin(theta.t.sample[i,],pois.data) - 
				dmt(theta.t.sample[i,], mean = c(tpar$m), S = tpar$var, df = tpar$df,log = TRUE)

wt.vec      	<- exp( d.vec - max(d.vec))
SIR.indices 	<- sample(1:(10^5),size=10^3,replace=T,prob=wt.vec)
theta.SIR.sample 	<- theta.t.sample[SIR.indices,]

plot(theta.SIR.sample[1:1000,1],theta.SIR.sample[1:1000,2])
mycontour(pois.loglin,c(1.5,4,-.2,.1),pois.data,add=T)

hist(theta.SIR.sample[,2], freq = FALSE) 
curve(dnorm(x,mean=fit$mode[2],sd = sqrt(fit$var[2,2])), 
				col = 2, lty = 2, lwd = 2, add = TRUE)

quantile(theta.SIR.sample[,2],prob=c(0.025,.5,.975))
fit$mode[2] ; 1.96 * sqrt(fit$var[2,2]) ; fit$mode[2] - 1.96 * sqrt(fit$var[2,2])


# Dalet 

plot(pois.data$i,pois.data$y,log="y")
lines(lowess(pois.data$i,pois.data$y,f = 1),col=2)

pred.laplace <- exp(fit$mode[1] + (1:18)*fit$mode[2])
lines(1:18,pred.laplace,col=3)					# What is this line???

lm.pois <- glm(pois.data$y ~ pois.data$i,family="poisson")
summary(lm.pois)

fit$mode
sqrt(diag(fit$var))



#	Q.2

# Bet

library(coda)


# (x,y) ~ N(theta.x,theta.y, (1,.8,.8,1))

#  x = 0, y = 0

x <- 0
y <- 0


theta.x <- rep(5, 100000)
theta.y <- rep(-5,100000)

for(i in 2:100000)
{
	theta.x[i] <- rnorm(1,mean = x + 0.8*(theta.y[i-1] - y),sd = sqrt(1 - 0.8^2))
	theta.y[i] <- rnorm(1,mean = y + 0.8*(theta.x[i]   - x),sd = sqrt(1 - 0.8^2))
}	

points(theta.x[1],theta.y[1],col=2,pch="+",cex=2)
lines(theta.x[1:50],theta.y[1:50],col=2)

t.test(theta.x[10001:20000])

t.test(theta.x[10001:20000])$conf[1:2]
t.test(theta.x[20001:30000])$conf[1:2]
t.test(theta.x[30001:40000])$conf[1:2]
t.test(theta.x[40001:50000])$conf[1:2]
t.test(theta.x[50001:60000])$conf[1:2]
t.test(theta.x[60001:70000])$conf[1:2]
t.test(theta.x[70001:80000])$conf[1:2]
t.test(theta.x[80001:90000])$conf[1:2]
t.test(theta.x[90001:100000])$conf[1:2]

theta.x	<- rnorm(100000)

t.test(theta.x[10001:20000])$conf[1:2]
t.test(theta.x[20001:30000])$conf[1:2]
t.test(theta.x[30001:40000])$conf[1:2]
t.test(theta.x[40001:50000])$conf[1:2]
t.test(theta.x[50001:60000])$conf[1:2]
t.test(theta.x[60001:70000])$conf[1:2]
t.test(theta.x[70001:80000])$conf[1:2]
t.test(theta.x[80001:90000])$conf[1:2]
t.test(theta.x[90001:100000])$conf[1:2]