Байесовский параметр c анализ выживаемости - PullRequest
Новая
       12

Байесовский параметр c анализ выживаемости

0 голосов
/ 18 марта 2020

Hello Stackoverflowers,

Я работал над уравнением, найденным в книге: Байесовский анализ выживаемости Джозефа Ибрагима 2001 (Глава параметри c модели, стр. 40-42). Мне удается заставить модель работать с усеченным гамма-распределением в R, но за всю жизнь я не понял, почему моя вероятность застрять близка к нулю. Ниже приведены мои коды для моделирования и выборки Гиббса, которые я кодировал.

Мое моделирование на основе параметризации пакета flexsurv: 

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)

Моя выборка Гиббса: 

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[,]

Любая помощь будет оценена.

Отчаянный программист,

...