In the following code I have created the posterior density for a Bayesian
survival model with four parameters. However, when I try to use the adapt
function to perform integration in four dimensions (on my old version of R
I get an error message saying that I have applied a non-function, although
the function does work when I type kernel2(param0, theta0), or on the
newer version of R the computer crashes.
I have attached the following R code:
hazard1 <- function(x, param) {
if(min(x, param)<=0) {return("function undefined")}
return(param[1]*x^(param[2] -1))
}
gm <- function( theta, param) {
dgamma(param[1], shape = .1, scale = .1) *
dgamma(param[2], shape = .1, scale = .1) *
dgamma(theta[1], shape = .1, scale = .1) *
dgamma(theta[2], shape = .1, scale = .1)
}
chazard1 <-function(x, param) {
if(min(x, param) <=0) {return("function undefined")}
param[1]*x^param[2] / param[2]
}
hazard2 <-function(x, theta, param) {
if(min(x, theta, param) <= 0) {return("function undefined")}
ratio = 1/(theta[1]*exp(-chazard1(x, param)) +
theta[2]*(1-exp(-chazard1(x, param))))
return(ratio*hazard1(x, param))
}
chazard2 <- function(x, theta, param) {
if(min(x, theta, param) <=0) {return("function undefined")}
val=c()
for (i in 1:length(x)) {
val=c(val, integrate(hazard2, lower=0, upper=x, theta=theta,
param=param)$value)}
return(val)
}
param0 = c(1.11, .833)
theta0 = c(1.11, .833)
yn1 <-read.table("yang1.txt", col.names = c("trtmt", "deathday", "reason"))
yn2 <-read.table("yang2.txt", col.names = c("trtmt", "deathday", "reason"))
for (i in 1:length(yn1$reason)) {
if (yn1$reason[i] == -1) (yn1$reason[i] = 0) else (yn1$reason[i] = 1) }
for (i in 1:length(yn2$reason)) {
if (yn2$reason[i] == -1) (yn2$reason[i] = 0) else (yn2$reason[i] = 1) }
x1 = yn1$deathday / (1*10^3);
d1 = yn1$reason;
x2 = yn2$deathday / (1* 10^3);
d2 = yn2$reason;
loglik <- function( theta, param) { sum(d1*log(hazard1(x1, param)) -
chazard1(x1,
param) + d2*log(hazard2(x2, theta, param)) - chazard2(x2, theta, param))}
/ 10000
kernel1 <- function(theta, param) { exp(loglik(theta, param)) + log(10000)}
kernel2 <- function(theta, param) {( exp(loglik(theta, param)) +
log(10000)) * gm(theta, param) }
kernel2(theta0, param0)
adapt(4, lower = c(0, 0, 0, 0), upper = c(1, 1, 1, 1), functn =
kernel2(theta0, param0))
So I am trying to integrate the posterior density with four different
parameters. But I cannot do anything to get the adapt function to work.
______________________________________________
R-help@stat.math.ethz.ch mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
