Bayesian parametric survival analysis

Question

Hello Stackoverflowers,

I have been working on the equation found in the book: Bayesian survival analysis by Joseph Ibrahim 2001 (Chapter parametric models p40-42). I manage to get a model going with a truncated gamma distribution in R but for the life of me, I have not figured out why my likelihood is stuck near zero. Below are my codes for both the simulation and the gibbs sampling that I coded.

My simulation based on flexsurv package parametrisation :

sim_gamma = function(N,shape,rate, val1, rateC,seed){
  set.seed(seed)
  X<-sort(rep(0:1,N/2))
  Time<-rgamma(N,shape=exp(shape), rate=exp((rate+val1*X))) 
  cens<-rgamma(N,shape=exp(shape), rate=exp(rateC)) 
  Y<-pmin(Time,cens)
  delta<-1*I(Time<cens)
  return(data.frame(time=Y,status = delta,x1 = X))
}

a = sim_gamma(N = 500,
              shape = 1.5, 
              rate = 0.3, 
              rateC =0.01,
              val1=0.5,
              seed =1)

My gibbs sampling :

library(cubature)
g <- function(u, alpha, lambda) {u^(alpha-1)*exp(-u)}
IG = function(yi, alpha, lambda){
  sapply(1:length(y), function(i)
    cubintegrate(g, lower = 0, upper = y[i]*exp(lambda), method = "pcubature", 
                 alpha = alpha, lambda= lambda)$integral)
}


y = a$time
status = a$status
d = sum(status)

# prior of alpha 
alpha0 =3
k0 = 0.01 #log(0.1) # same format as lambda

# prior of lambda 
sigma0 = 1
mu0 = 0


alpha = exp(1.5) #shape > 0
rate = exp(0.1) # rate > 0
lambda = log(rate)

################## ## ## ## ## ## ## ## ## ## ## ## ## ##################

n = 20
dat = matrix(NA, nrow=n, ncol=4)

alpha = 2 #shape > 0
rate = 0.5 # rate > 0
lambda = log(rate)
y_star = y
y_star = (y_star^(alpha-1)) / (gamma(alpha)*((1-1/gamma(alpha)*IG(y[i],alpha,lambda)))) * exp(alpha*lambda - y_star * exp(lambda))
y_star[status==1] = y[status==1] 
event = length(y)


for(i in 1:n){
  Lik = ((1/gamma(alpha))^event) *exp(alpha*lambda*event + sum(y_star*(alpha-1)  - y_star*exp(lambda))) 
  joint_posterior = Lik * alpha^(alpha0-1) *exp(- k0*alpha - 1/(2 *sigma0^2)*(lambda - mu0)^2 )
  alpha = Lik * alpha^(alpha0-1)*exp(-alpha0*k0)
  lambda = Lik *exp(-1/2*((lambda - mu0) / sigma0)^2)
  dat[i,]=(c(Lik,alpha, lambda,exp(lambda)) )
  y_star = (y_star^(alpha-1)) / (gamma(alpha)*((1- 1/gamma(alpha)*IG(y[i],alpha,lambda)))) * exp(alpha*lambda - y_star * exp(lambda))
}
dat[,]

Any help would be appreciated.

A desperate programmer,