# Code for problem 3 of hw 3

y<-c(0,1,3,5)
n<-c(5,5,5,5)
dose<-c(-0.863,-0.296,-0.053,0.727)
mu <- c(0.8732163,7.9110577)
Sigma <- matrix(data=c(1.032047^2,0.7251766*1.032047*4.98319,
0.7251766*1.032047*4.98319, 4.98319^2),nrow=2,ncol=2)

# Sample from the approximating Normal distribution and scatterplot

par(mfrow=c(2,2))
m <- 2000
alphaN <- c(m,0,0)
betaN <- c(m,0,0)
for (i in 1:m) {
x <- mu + t(chol(Sigma)) %*% rnorm(2)
alphaN[i] <- x[1]
betaN[i] <- x[2]
}
plot(alphaN,betaN)
LD50 <- -alphaN/betaN
LD50 <- LD50[LD50>-1.0 & LD50 <1.0] # to avoid extreme sample values of LD50
                                    # that distort the histogram
hist(LD50, xlab="sample from appr. posterior density ofLD50",
nclass=30,prob=TRUE)

#importance resampling (SIR)

posterior<-function(a,b){(exp(a +b*(-0.863))/(1+exp(a+ b*(-0.863))))^0*
(1- exp(a +b*(-0.863))/(1+exp(a +b*(-0.863))) )^5 *
(exp(a +b*(-0.296))/(1+exp(a +b*(-0.296))))^1*
(1- exp(a +b*(-0.296))/(1+exp(a+b*(-0.296))) )^4 *
(exp(a +b*(-0.053))/(1+exp(a+b*(-0.053))))^3*
(1- exp(a+b*(-0.053))/(1+exp(a+b*(-0.053))) )^2 *
(exp(a +b*(0.727))/(1+exp(a+b*(0.727))))^5*
(1- exp(a +b*(0.727))/(1+exp(a+b*(0.727))) )^0
}

proposal<-function(a,m,S){
1/(2*pi*sqrt(det(S)))*exp(-1/2*t(a-m)%*%solve(S)%*%(a-m))
}
w<-seq(m)
for (i in 1:m) {
ab<- c(alphaN[i],betaN[i])
w[i] <- posterior(alphaN[i],betaN[i])/proposal(ab,mu,Sigma)
}
w <- w/sum(w)

hist(w,xlab="sample of importance weights", nclass=50,prob=TRUE)

index <- sample(seq(m),size=500,replace=T,prob=w)
for (i in 1:m){
alphaSIR[i] <- alphaN[index[i]]
betaSIR[i] <- betaN[index[i]]
}
plot(alphaSIR,betaSIR)