# 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);

//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."

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

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