Defining Jeffreys' prior for transformed parameters and censored data

Good day all,

I am writing a code for Bayesian analysis with transformed parameters of the Weibull distribution. I have to define a Jeffreys’ prior to one of the transformed variables (\rho) and the transformations are non-linear. I have tried working around the Jacobian adjustment for the transformed parameters. Could anyone please help me with identifying any mistakes in the model? I am supposed to use this coding for a simulation study with censored data.

Weibull Parameters
\alpha - scale
\beta - shape

transformed parameters
\mu = log(\alpha)
\rho = 1/(\beta*log(\alpha))

Jeffreys’ prior for rho is 1/\rho

My Stan model is as follows.

model_Jeffreys <- "


data {
int<lower=0> J; // number of observations
int<lower=0> N_cen; // number of censored observations
vector [J]sampledata; // failure in hours
vector [N_cen]y_cen; //values of censored data


}

parameters {
real alpha; // scale
real<lower=0> beta; // shape


}


transformed parameters{
real mu;
real<lower = 0> rho;

rho <- 1 /(beta*log(alpha));

mu <- log (alpha);


}

model {
//prior
mu ~ uniform (0, 20);
increment_log_prob(-log(mu)); //adjustment for Jacobian


//Likelihood and prior
sampledata ~ weibull (1/(rho*mu), exp(mu));
increment_log_prob(weibull_ccdf_log (y_cen,  1/(rho*mu), exp(mu)));
increment_log_prob(-log(rho));


}
"
mod_Jeffreys <- stan_model(model_code = model_Jeffreys)

Thank you.

Hi!
maybe try to use this:

data {
   int<lower=0> J;            // number of observations
   int<lower=0> N_cen;   // number of censored observations
   vector[J] sampledata;  // failure in hours
   vector[N_cen] y_cen;  // values of censored data
}
parameters {
   real mu;                       // scale
   real<lower=0> rho;      // shape
}
model {
  //prior
  target += uniform_lpdf(mu|0,20);
  target += -log(rho);  
  //Likelihood
  target += weibull_lpdf(sampledata | 1/(rho * mu), exp(mu));
  target += weibull_lcdf(y_cen |  1/(rho * mu), exp(mu));
}

I haven’t run it myself but I worked with a Jeffreys prior some days ago and I solve it similarly. If you still have some trouble I will help and debug myself the stan code tomorrow morning (I have homework now) :)

Stan has some recommendations for prior selection, because of Jeffrey’s prior work really bad, if you haven’t read it check here, could be helpful.

2 Likes

Thank you for the reply. I tried using the above code and got the following error. Earlier when I tried to use target +=, I got the same error and thus, proceeded with increment\_log\_prob.

"SYNTAX ERROR, MESSAGE(S) FROM PARSER:

variable “target” does not exist.

ERROR at line 14

12: model {
13: //prior
14: target += uniform_lpdf(mu|0,20);
^
15: target+= -log(rho);

Error in stanc(file = file, model_code = model_code, model_name = model_name, :
failed to parse Stan model ‘c5e7c7c843af9f2827edfff1ef4dc754’ due to the above error."

Could you please help me resolving this? Thank you.

Do tou have the data to try it?

Here’s the code for data generation and a test data set.

#Simulate data
weibull.data <- function (samplesize, #sample size
                          scale, #scale parameter
                          rho,
                          propFailed){
  data.wei <- sort (rweibull (samplesize, shape = 1/(scale*rho), 
                              scale = exp(scale)))
  T.cen <- data.wei[max(2, ceiling(propFailed*samplesize))]
  censored <- data.wei > T.cen
  event <- c (rep(1, length(data.wei[censored == FALSE])), rep(0, length(data.wei[censored == TRUE])))
  gen.fails <- c (data.wei[censored == FALSE], rep(NA, length(data.wei[censored == TRUE])))
  rank.fails <- rank (gen.fails, ties.method = "first", na.last = "keep")
  gen.fails [is.na(gen.fails)] <- T.cen
  y.cen <- rep (T.cen, length(data.wei[censored == TRUE]))
  
  
  simul.data <- data.frame(sampledata = data.wei,
                           Censored = censored,
                           obs.fails = gen.fails,
                           ranks = rank.fails)
  
  sim.data.bay <- list (sampledata = data.wei[censored == FALSE],
                        obs_fails = gen.fails,
                        Censoreddata = censored,
                        Censored = event,
                        J = length (data.wei[censored == FALSE]),
                        N_cen = length(data.wei[censored == TRUE]),
                        y_cen = y.cen)
  
  return(sim.data.bay)
}



set.seed(1589)
(test <- (weibull.data(samplesize = 50,
                      scale = 1,
                      rho = 0.25,
                      propFailed = 0.25)))

Thank you very much for your help @asael_am

Hello, was a simple syntax error, just in case I attach the stan file here: stantry.stan (521 Bytes)

library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
Sys.setenv(LOCAL_CPPFLAGS = '-march=native')

set.seed(1589)
test = (weibull.data(samplesize = 50,scale = 1,rho = 0.25,propFailed = 0.25))
smodel = stan_model(file = "stantry.stan")
sf1 = sampling(object = smodel,data = test,chains = 1,iter = 2000)
sf1

It runs ok, and the code above is also corrected, so I hope it works for you :)

1 Like

I still get the same error with the variable target.

what stan version are you using? I am using rstan 2.19.3 and it works ok this is what sf1 throws to me:
> sf1

Inference for Stan model: stantry.
1 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=1000.

       mean se_mean   sd   2.5%    25%    50%   75% 97.5% n_eff Rhat
mu     0.61    0.00 0.04   0.51   0.58   0.61  0.63  0.67   365    1
rho    0.26    0.00 0.08   0.15   0.20   0.25  0.30  0.45   332    1
lp__ -10.37    0.06 1.10 -13.05 -10.73 -10.05 -9.58 -9.32   385    1

Samples were drawn using NUTS(diag_e) at Thu Apr 02 09:26:03 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

The only thing I need to fix was this:

For this:

target += -log(rho);

And even so I dont think this is the real problem :(

1 Like

I am using rstan 2.9. May be that is the problem. I will try updating the package :).

do you manage to solve it?

Yes @asael_am. I used increment_log_prob instead of target as I was getting the same error. Thank you for your help.

1 Like

Ok nice! Glad it works for you!

1 Like