# Same data as before

Y <- c(65, 156,100,134,16,108,121,4,39,143,56,26,22,1,1,5,65)
wbc <- c(2300,750,4300,2600,6000,10500,10000,17000,5400,7000,9400,
	 32000,35000,100000,100000,52000,100000)

x <- log(wbc/10000)

# logposterior must be a function of theta only; other datavalues are
# specified using the data option.

logpost <- function(theta) {
  theta0 <- theta[1]
  theta1 <- theta[2]
  mu <- theta0*exp(- theta1*x)
  ll <- sum(log(dexp(y, 1/mu)) )
return(ll) }

library(sbayes)   #  currently only works on the DEC stations ?
                  #  sbayes is available from statlib.
                  #  http://lib.stat.cmu.edu/

# Usage:  function for logposterior, starting values, scale, other data
lp <- bayes.model(logpost,c(35,.8),scale=c(10, .5), data=list(x=x,y=Y))

bayes.moments(lp,method="1st")

# Laplace or second order approximation
bayes.moments(lp, method="2nd")

# option covar=T will give 2nd order covariance matrix of theta
bayes.moments(lp, method="2nd", covar=T)

motif()
plot(bayes.mar1(lp, 1, seq(20,130,len=15)), 
	type="l", xlab="theta0", ylab="P(\theta0 | Y)") 

lines(seq(20,130,len=50), dnorm(seq(20,130,len=50), mean=56.8, sd=13.97), lty=2)

# P(Y* > 52 | Y) = probability of survivng more than 1 year if wbc=50,000

g <- function(theta) {
	mu <- theta[1]* exp(-log(50000/10000) * theta[2])
	exp(-52 /mu) }

bayes.moments(lp,g, method="1st")    

bayes.moments(lp,g, method="2nd")                        

plot(bayes.mar1(lp, g, seq(.00,.3, len=10), method="2nd"), type="l", 
	xlab="1 Year Survival Probability",
	ylab="Posterior Density")

lines(seq(0,.3, len=100), dnorm(seq(0,.3, len=100), mean=.18, sd=.11), lty=2)