N <- 10000
burnin <- 1000
k <- 10
y <- c(5,1,5,14,3,19,1,1,4,22)
t <- c(94.320,15.720,62.880,125.760, 5.240,
       31.440, 1.048, 1.048,  2.096,10.480)
r <- y/t
U <- 1.0
rho <- 0.2
nu <- 1.4
S <- rep(NA,N)
S[1] <- 2.0
theta <- matrix(NA,nrow=N,ncol=k)
theta[1,] <- rep(1.0,k)
for (j in 2:N) {
  for (i in 1:k) {
    theta[j,i] <- (S[j-1]+2*t[i])^(-1)*rchisq(1,nu+2*y[i])
  }
  S[j] <- (U+sum(theta[j,]))^(-1)*rchisq(1,rho+k*nu)
}
Strunc <- S[burnin:N]
thetatrunc <- theta[burnin:N,]
thetamean <- apply(thetatrunc,2,mean)
thetasd <- apply(thetatrunc,2,sd)
thetastats <- cbind(thetamean,thetasd)
colnames(thetastats) <- c("mean","sd")
cat("\nPump means estimated as\n")
print(thetastats)
cat("\nS has mean",mean(Strunc),
    "and s.d.",sd(Strunc),"\n")