Unable to sample from the prior distribution

Modeling
1

Hi,
I’m working with the following model. I tried to run it and it works. Now I want to sample from the prior distribution by commenting out sor ~ beta_proportion(inv_logit(mu), kappa) . See below.

data {
  int<lower=1> N_obs;
 
  vector[N_obs] sor;
  vector[N_obs] net_dist_km;
  vector[N_obs] pp_diff;
  
  int<lower=1> N_flow;
  vector[N_obs] flow_conn;
 
  int<lower=1> N_strahlers;
  vector[2] alpha;
  int<lower=1, upper=N_strahlers> strahler_diff[N_obs];
  
  int N_basins;
  int<lower=1, upper=N_basins> basin_idx[N_obs];
 
  int N_samples;
  int<lower=1, upper=N_samples> sample_idx1[N_obs];
  int<lower=1, upper=N_samples> sample_idx2[N_obs];

}

transformed data {
  vector[N_obs] delta_net_dist_km = (net_dist_km - 100.0) / 100.0;
  vector[N_obs] delta_pp_diff = (pp_diff - 5.71) / 5.71;
}

parameters {

  real a_baseline;
  
  vector[N_basins] a_basin;
  
  vector[N_samples] a_sample;

  vector[N_basins] bstrah; 
  real mu_bstrah;
  real<lower=0> sigma_bstrah;
  simplex[2] delta;


  vector[N_basins] bflow;
  real mu_bflow;
  real<lower=0> sigma_bflow;
  
  vector[N_basins] bnet;
  real mu_bnet;
  real<lower=0> sigma_bnet;

  vector[N_basins] bpp;
  real mu_bpp;
  real<lower=0> sigma_bpp;
  
  real<lower=0> kappa;
}  

transformed parameters {
  vector[N_obs] mu; 
  vector[N_obs] delta_sum; 
  vector[3] delta_j = append_row(0, delta);
  
  
   for ( i in 1:N_obs ) {
   delta_sum[i] = sum(delta_j[1:strahler_diff[i]]);
   }

 
  mu =  a_baseline
      + a_sample[sample_idx1] + a_sample[sample_idx2]
      + a_basin[basin_idx] 
      + bstrah[basin_idx].*  delta_sum
      + bflow[basin_idx] .* flow_conn
      + bnet[basin_idx] .* delta_net_dist_km
      + bpp[basin_idx] .* delta_pp_diff;
}

model {
  a_baseline ~ normal(0, 1);
  a_sample ~ normal(0, 1);

  bstrah~normal(mu_bstrah,sigma_bstrah);
  mu_bstrah~normal(0,1);
  sigma_bstrah ~ exponential(1);
  delta ~ dirichlet(alpha);
  

  bflow~normal(mu_bflow,sigma_bflow);
  mu_bflow~normal(0,1);
  sigma_bflow ~ exponential(1);
  
  bnet ~ normal(mu_bnet, sigma_bnet);
  mu_bnet~normal(0,1);
  sigma_bnet~exponential(1);
  
  bpp ~ normal(mu_bpp, sigma_bpp);
  mu_bpp ~ normal(0,1);
  sigma_bpp ~ exponential(1);
  
  kappa ~ normal(0, 50);

  //sor ~ beta_proportion(inv_logit(mu), kappa);
}

generated quantities{
  real sor_prior[N_obs] = beta_proportion_rng(inv_logit(mu), kappa);
}

However, STAN throws this error:

recompiling to avoid crashing R session
starting worker pid=10728 on localhost:11073 at 12:07:09.789
starting worker pid=10742 on localhost:11073 at 12:07:09.947
starting worker pid=10756 on localhost:11073 at 12:07:10.105
starting worker pid=10770 on localhost:11073 at 12:07:10.270

SAMPLING FOR MODEL 'ppd' NOW (CHAIN 1).

SAMPLING FOR MODEL 'ppd' NOW (CHAIN 2).
Chain 1: 
Chain 1: Gradient evaluation took 0.001148 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 11.48 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'ppd' NOW (CHAIN 3).
Chain 2: 
Chain 2: Gradient evaluation took 0.002042 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 20.42 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)

SAMPLING FOR MODEL 'ppd' NOW (CHAIN 4).
Chain 3: 
Chain 3: Gradient evaluation took 0.001202 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 12.02 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 4: 
Chain 4: Gradient evaluation took 0.001402 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 14.02 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 3: Exception: beta_rng: Second shape parameter is 0, but must be > 0!  (in 'model765e309fc911_ppd' at line 106)

Chain 3: Exception: beta_rng: Second shape parameter is 0, but must be > 0!  (in 'model765e309fc911_ppd' at line 106)

Chain 3: Exception: beta_rng: Second shape parameter is 0, but must be > 0!  (in 'model765e309fc911_ppd' at line 106)

Chain 3: Exception: beta_rng: Second shape parameter is 0, but must be > 0!  (in 'model765e309fc911_ppd' at line 106)

Chain 3: Exception: beta_rng: Second shape parameter is 0, but must be > 0!  (in 'model765e309fc911_ppd' at line 106)

Chain 3: Exception: beta_rng: Second shape parameter is 0, but must be > 0!  (in 'model765e309fc911_ppd' at line 106)

Line 106 corresponds to the generated quantities section. What am I doing wrong?

Thanks

2

Probably you did nothing wrong. This happens when inv_logit(mu) or 1 - inv_logit(mu) is close to 0.
You can try
beta_rng(inv_logit(mu) * kappa, kappa - kappa * inv_logit(mu))
or try
beta_rng(exp(log_inv_logit(mu) + log(kappa)), exp(log1m_inv_logit(mu) + log(kappa)))

Mostly it happens in warm-up phase, you might write an exception handling

if(parameter2 == 0)
sampling = 1;

3

Thanks

I added:

generated quantities{
real sor_prior[N_obs] = beta_proportion_rng(inv_logit(mu) * kappa, kappa - kappa * inv_logit(mu));
}

but got:

Chain 4: Exception: beta_proportion_rng: Location parameter[8] is 1.91111, but must be less than or equal to 1  (in 'model765e788f8645_ppd' at line 106)

Chain 4: Exception: beta_proportion_rng: Location parameter[8] is 1.91111, but must be less than or equal to 1  (in 'model765e788f8645_ppd' at line 106)

Chain 4: Exception: beta_proportion_rng: Location parameter[8] is 1.91111, but must be less than or equal to 1  (in 'model765e788f8645_ppd' at line 106)

Chain 4: Exception: beta_proportion_rng: Location parameter[8] is 1.91111, but must be less than or equal to 1  (in 'model765e788f8645_ppd' at line 106)

Did you notice I’m using beta_proportion? Maybe it requires a different formulation?

Concerning exceptions: I’m still new to STAN. Where should I add such an exception? And should it be:

if(mu == 0)
sampling = 1;

?

4

beta_rng not beta_proportion_rng.

beta(a, b), just check if one of the parameters is 0 and set the sampling parameter accordingly.

5

real sor_prior[N_obs] = beta_rng(inv_logit(mu) * kappa, kappa - kappa * inv_logit(mu));

causes:

Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
Chain 2: Exception: beta_rng: Second shape parameter[3421] is 0, but must be > 0!  (in 'model765e3417a1f9_ppd' at line 106)

Chain 2: Exception: beta_rng: Second shape parameter[3421] is 0, but must be > 0!  (in 'model765e3417a1f9_ppd' at line 106)

Chain 1: Exception: beta_rng: Second shape parameter[57] is 0, but must be > 0!  (in 'model765e3417a1f9_ppd' at line 106)
6

Same? I don’t exactly know where inv_logit equals to 1.

1/(1+exp(-16))
[1] 0.9999999

if(mu>16) sampling = 0;
if(mu<-16) sampling = 1;

7

Yes, same problem:

Chain 1: Exception: beta_rng: First shape parameter[1755] is 0, but must be > 0!  (in 'model765e5369892a_ppd' at line 106)

Chain 1: Exception: beta_rng: First shape parameter[1755] is 0, but must be > 0!  (in 'model765e5369892a_ppd' at line 106)

Chain 1: Exception: beta_rng: First shape parameter[1755] is 0, but must be > 0!  (in 'model765e5369892a_ppd' at line 106)

Sorry, but where exactly should I add the exception? Example of code that does not work:

    generated quantities{
    if(mu>16) sampling = 0;
    if(mu<16) sampling = 1;
    real sor_prior[N_obs] = beta_proportion_rng(inv_logit(mu),kappa);
    }
8 
generated quantities{
real sor_prior[N_obs];
for(i in 1:N_obs) {
if(mu[i] >16) sor_prior[i] = 0;
else if(mu[i] < -16) sor_prior[i] = 1;
else
sor_prior[i] = beta_proportion_rng(inv_logit(mu[i]),kappa);
}}
9

Thank you.

1 Like