#source("f:\\David\\Enseignements\\Cours ICASA\\MiniBayes\\ISgrowth.txt")

####Estimation of growth rate using an Importance Sampling algorithm####

###Seed used for random sampling###
set.seed(1)

###Number of simulations###
N<-10000

###Prior distribution parameters###
mu<-0.03
tau<-0.015

###Data###
Y<-c(0.83,0.95)
X<-c(100,200)
sigma<-0.02

###Simulations from prior###
thetaPrior<-rnorm(N,mu,tau)
Likeli<-matrix(nrow=length(thetaPrior),ncol=length(Y))
LikeliProd<-thetaPrior
#par(mfrow=c(1,1))
#plot(thetaPrior, dnorm(thetaPrior, mu, tau), cex=1, xlab="Theta", ylab="Density")

###Weight calculation###
for (i in 1:N) {
	for(j in 1:length(Y)) {
			Model.i.j<-1-exp(-thetaPrior[i]*X[j])
			Likeli[i,j]<-dnorm(Y[j],Model.i.j,sigma)
							}
			LikeliProd[i]<-prod(Likeli[i,])				
					}
					
Weight<-LikeliProd	

###Weight normalization###
Weight<-LikeliProd/sum(LikeliProd)

###Plot of posterior###
par(mfrow=c(2,2))
plot(thetaPrior,Weight,xlab="Theta",ylab="Weight", type="l")
title("Plot of normalized weights")

###Resampling###
thetaPost<-sample(thetaPrior,replace=T,prob=Weight)
hist(thetaPost,labels=F,xlab="Theta",main="Parameter values drawn from posterior")	

print("Posterior mean:")	
print(mean(thetaPost))
print("Posterior variance:")
print(var(thetaPost))