Defining Jeffreys' prior for transformed parameters and censored data

Modeling
,
1

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.

2

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
3

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.

4

Do tou have the data to try it?

5

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

6

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
7

I still get the same error with the variable target.

8

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).
9

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
10

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

11

do you manage to solve it?

12

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
13

Ok nice! Glad it works for you!

1 Like