##Example of importance sampling method
##Exercise 10.9, and 10.10 pag 319, Gelman's book
importance(w=1,10000,1000,mmm=0,sss=2,DF=3,xl=c(-6,6))
importance2(w=1,10000,1000,mmm=0,sss=2,DF=3,xl=c(-6,6))



importance_function(w=0,L,K,mmm=0,sss=1,DF=3,xl=c(-4,4)){
  if (w == 1) postscript("/home/biostats/fdominic/course/SIRgood.ps")
  par(mfrow=c(2,1))
  par(oma=c(0,0,2,0))
  theta        _  NULL
  theta.better _ NULL  
  weight       _ NULL 
  for(l in 1:L){ 
    theta[l]    _ sss*rt(1,DF)
    weight[l]   _ dnorm(theta[l],mean=mmm,sd=sss)/dt(theta[l],df=DF)
  }
  hist(log(weight),nclass=100,xlab="log importance ratios",density=-1)
  post.mean    _ mean(theta*weight)
  post.var     _ 1/L*sum(weight*(theta-mean(theta))^2)
  theta.better _ sample(theta, K,prob=weight)
   hist(theta.better,xlab="theta",nclass=30,
       xlim=c(xl[1],xl[2]) ,density=-1,prob=T,yaxt="n")
  abline(v=mean(theta.better))
  abline(v=mean(theta.better)-sqrt(var(theta.better)),lty=2)
  abline(v=mean(theta.better)+sqrt(var(theta.better)),lty=2)
  par(mfrow=c(1,1))
  par(oma=c(0,0,0,0)) 
  if (w == 1) dev.off()
  return(post.mean,post.var)
}


importance2_function(w=0,L,K,mmm=0,sss=2,DF=3,xl=c(-5,5)){
  if (w == 1) postscript("/home/biostats/fdominic/course/SIRbad.ps")
  par(mfrow=c(2,1))
  par(oma=c(0,0,2,0))
  theta_NULL  
  weight_NULL 
  theta.better_NULL
  weight_NULL 
  for(l in 1:L){ 
  theta[l]  _ rnorm(1,mean=mmm,sd=sss)
  weight[l]  _ dt(theta[l],df=DF)/dnorm(theta[l],mean=mmm,sd=sss)}
  hist(log(weight),nclass=100,xlab="log importance ratios",density=-1)
  post.mean _ mean(theta*weight)
  post.var  _ 1/L*sum(weight*(theta-mean(theta))^2)
  theta.better _ sample(theta, K,prob=weight)
  hist(theta.better,xlab="theta",nclass=100,
       xlim=c(xl[1],xl[2]) ,density=-1,prob=T,yaxt="n")
  abline(v=mean(theta.better))
  abline(v=mean(theta.better)-sqrt(var(theta.better)),lty=2)
  abline(v=mean(theta.better)+sqrt(var(theta.better)),lty=2)
  par(mfrow=c(1,1))
  par(oma=c(0,0,0,0)) 
  if (w == 1) dev.off()
 return(post.mean,post.var)
}