Generating random numbers from a conditional distribution

Modeling
1

Hello,

Everyone! I am new to Stan and am trying to use it to fit a Weibull distribution. I got my model correct to get posterior estimates. Now I need to get some predicted values based on a conditional Weibull, i.e., conditional on the predicted value is greater than say ti (i=1, …, q_n) for q_n subjects. I tried the generated quantities block below, but it takes forever to run, how can I simplify it?

data {
  int<lower=0> n;                                // number of subjects
  real adt[n];                                        // survival time
  int<lower=0, upper=1> evt[n];          // event indicator
  int<lower=0> q_n;                            // number of subjects with administrative censor
}
parameters {
  real<lower=0> alpha;          
  real<lower=0> sigma;    
}
model {
  for (i in 1:n){
    target += weibull_lpdf(adt[i]|alpha, sigma)*evt[i]+weibull_lccdf(adt[i]|alpha, sigma)*(1-evt[i]); 
  }
}
generated quantities{
  vector[q_n] adt_pred;
  for (j in 1:q_n){
    adt_pred[j] =0;
    while(adt_pred[j]<adt_q[j]){
      adt_pred[j] = weibull_rng(alpha, sigma);
      }
  }
}
1 Like
2

hello again,

I am able to create a function to make generate the conditional weibull, with the code below:

functions {
  real weibull_trunc_rng(real alpha, real sigma, real lb){
    real lp = weibull_cdf(lb, alpha, sigma);
    // real up = weibull_cdf(ub, alpha, sigma);
    real u = uniform_rng(lp, 1);
    real y = sigma * (-log1m(u))^inv(alpha);
    return y;
  }
}
...
generated quantities{
  vector[q_n] adt_pred;
  // vector[q_n] adt_pred_order;
  // real predt;
  for (j in 1:q_n){
    adt_pred[j] = weibull_trunc_rng(alpha, sigma, adt_q[j]);
  }
  // adt_pred_order= sort_asc(adt_pred);
  // predt=sort_asc(adt_pred)[addevt];
}

But when I add the commented out codes to generate more, i got the following error message:
Error in sampler$call_sampler(args_list[[i]]) :
Exception: Exception: uniform_rng: Upper bound parameter is 1, but must be greater than 1 (in ‘model3e5f28701a7_evt_pred_weibull’ at line 5)
(in ‘model3e5f28701a7_evt_pred_weibull’ at line 33)

Any help on what is causing this and how to solve it will be appreciated!!

1 Like
3

Hi @Lei_Hua, welcome to the Stan forums! Would it be possible to share the data you’re using (or some simulated data that also works) as well as a the full Stan program that results in the error? I think that will help me or someone else debug this. Thanks!

4

Sure, I can provide the data. What is the easiest way of providing the data?

This happens when I work on something similar. I think the issue is about the function generating u.

5

Thanks. When responding to this post, one of the editing-related buttons at the top will let you upload a csv file. Or you can just share some R or Python code for generating data that demonstrates the same problem. Whatever you prefer.

6

Also one thing I just noticed is that adt_q doesn’t seem to be defined anywhere in the code you shared (unless I missed it). I see adt and adt_pred, but there’s an adt_q variabled used that I can’t find the definition of.

7

Thanks, @jonah, I have the adt_q defined the data block. I find a misspecification of my model. Fixing the models seems solved the problem now.

1 Like
8

Ok great, glad you solved it!

1 Like
9

My issue comes back, i’ve uploaded the data as an R objectmydata.csv (24.4 KB)

functions {
  real exponential_trunc_rng(real lambda, real lb){// function to generate truncated exponential;
    real lp = exponential_cdf(lb, lambda);
    real u = uniform_rng(lp, 1);
    real y = -log1m(u)/lambda;
    return y;
    }
  int fst_index(vector x, real y) {// function to find the first positio of matches;
    real tmp = x[1];
    int pos = 1;
    while(tmp<y){// to find location
      pos += 1;
      tmp = x[pos];
      }
    return pos;
    }
  }
data {
  int<lower=0> n;
  real adt[n];
  int<lower=-1, upper=1> evt[n];
  int<lower=0> q_n;
  real<lower=0> adt_q[q_n];
  real<lower=0> enrdt[q_n];
  real<lower=0> addevt;
}
parameters {
  // for weibull event time;
  real<lower=0> lambda_e;
  // for exponential discontinuation time;
  real<lower=0> lambda_d;
}
model {
  for (i in 1:n){
    target += exponential_lpdf(adt[i]|lambda_e)^(evt[i]==1) + exponential_lccdf(adt[i]|lambda_e)^(evt[i]!=1) +
              exponential_lpdf(adt[i]|lambda_d)^(evt[i]==-1) + exponential_lccdf(adt[i]|lambda_d)^(evt[i]!=-1);  
              }
}
generated quantities{
  real evt_pred[q_n];
  real disc_pred[q_n];
  real adt_pred[q_n];
  real evtdt_pred[q_n];
  real evtdt_pred_order[q_n];
  vector[q_n] ind;
  vector[q_n] ind_order;
  real predt;
  for (j in 1:q_n){
    evt_pred[j] = exponential_trunc_rng(lambda_e, adt_q[j]);
    disc_pred[j] = exponential_trunc_rng(lambda_d, adt_q[j]);
    adt_pred[j] = fmin(evt_pred[j], disc_pred[j]);
    ind[j] = (evt_pred[j]<=disc_pred[j]);
    evtdt_pred[j] = enrdt[j] + adt_pred[j];
  }
  evtdt_pred_order=evtdt_pred[sort_indices_asc(evtdt_pred)];
  ind_order=ind[sort_indices_asc(evtdt_pred)];
  // predt = evtdt_pred_order[min(which_equal(cumulative_sum(ind_order), addevt))];
  predt = evtdt_pred_order[fst_index(cumulative_sum(ind_order), addevt)];
}

below is my R code to get this work in R

mydata_mod <- list(
n = nrow(mydata),
adt = mydata$adt,
evt = mydata$evt,
# to simulate posterior dist and get estimates;
q_n = sum(mydata$evt==0),
adt_q = mydata$adt[mydata$evt==0],
enrdt = mydata$enrdt[mydata$evt==0],
addevt = 118-100
)
fit.exp <- stan(
file = “/usrfiles/shared-projects/clinicalInformatics/rstan/LeiHua_demo/Lei’s code/event prediction/bayesian prediction/oneArm/draft/evt_pred_exp.stan”, # Stan program
data = mydata_mod, # named list of data
# seed=1234,
chains = 4, # number of Markov chains
warmup = 1000, # number of warmup iterations per chain
iter = 5000, # total number of iterations per chain
cores = 2, # number of cores (could use one per chain)
refresh = 0 # no progress shown
)

When I run this with 4 chains, 3 chains do not have samples, while 1 has. It gave me the following errors:

[1] “Error in sampler$call_sampler(args_list[[i]]) : "
[2] " Exception: Exception: uniform_rng: Upper bound parameter is 1, but must be greater than 1 (in ‘model68ff229467f6_evt_pred_exp’ at line 4)”
[3] " (in ‘model68ff229467f6_evt_pred_exp’ at line 50)"

It seems the function I defined using runif is causing such. How can I modify to make it work?

10

Sorry, I just saw this. I think this was solved in Error with uniform_rng, right? If so, that’s great. If not then let me know and I can take another look.

11

Yes, thanks!