R help - Nov 2006 - Bayesian question (problem using adapt)

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.

On Sat, Nov 11, 2006 at 11:19:46AM -0500, cmdicaso at ncsu.edu wrote:
>  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))
I could not walk through your code because you failed to create a
simple reproducible example (eg the data is missing), but
kernel2(theta0,param0) is no longer a function --- kernel2 is, you
need to just supply that.

A couple of remarks:

- instead of "function undefined", you might just want to return NA or
NaN as is standard in R.

- why do you multiply by 10000?

- are you just integrating the posterior?  why?

- have you considered MCMC?  On state spaces bigger than 2 dimension,
  it is usually preferable to numerical integration (but since I don't
  know your model in detail, there might be something special about
  it).

Best,

Tamas

Is it possible there is a typo in your attached code?  You wrote... 

loglik <- function(theta, param) {
  sum(d1*log(hazard1(x1, param)) - chazard1(x1,param) + 
  d2*log(hazard2(x2, theta, param)) - chazard2(x2, theta, param))
} / 10000


Confidentiality Notice: This e-mail message, including any a...{{dropped}}