# R code for generating samples from the joint posterior density
# of (mu,sigma2) in Normal(mu,sigma2) data model with non informative
# prior p(mu,sigma2) proportional to 1/sigma2.
#
norm2.r_function(m)
{
 data_read.table("norm2.dat",header=T)
 surv_data[,1]
 y_log(surv)
 n_length(y)
 # s2: sample variance
 s2_1/(n-1)*sum((y-mean(y))^2)

# Simulation from the marginal posterior density for sigma2 
 chi_rchisq(m,n-1)
 sigma.post_(n-1)*s2/chi

# Simulation from conditional posterior density for mu 
# and from the (conditional) posterior predictive density 
# for y.new. 
# Simulation from a transformation of y.new: the probability that 
# survival time exceedes 150 (log150=5.01) 
 
 mu.post.cond_c()
 y.new_c()
 y.bar_mean(y)
 g.150_c()
 for(i in 1:m){
 mu.post.cond_c(mu.post.cond,rnorm(1,y.bar,sqrt(sigma.post[i]/n)))
 y.new_c(y.new,rnorm(1,mu.post.cond[i],sqrt(sigma.post[i]))) 
 g.150_c(g.150,as.numeric(y.new[i]>5.01))}

# Alternative: simulation from posterior predictive distribution for 
# a future observation y.new; probability computation of g.150 mean 
# value
#
# ti_rt(10000,19)
# y.new_sqrt(1+1/n)*sqrt(s2)*ti+y.bar 
# Posterior probability of survival time exceeding 150 
# g.150_1-pt((log(150)-y.bar)/(sqrt(1+1/n)*sqrt(s2)),19)

return(list(sigma.post,mu.post.cond,y.new,g.150))
}