Truncating a cauchy_rng() sampler

Modeling
1

I am trying to sample from the prior distribution in rstan using the generated quantities block (see https://stackoverflow.com/questions/57703920/sampling-from-prior-without-running-a-separate-model/57718910#57718910). For the prior on the noise distribution I wanted to use a cauchy distribution truncated at zero so it only generates positive values. In the accepted answer to my original post @Bob Carpenter suggested using the cauchy_rng() function within the generated quantities block.

Unfortunately however truncating the cauchy_rng() monte carlo generator results in an error. The model works if you don’t truncate the cauchy_rng() sampler but it generates values below 0.

Here’s toy data and the model to illustrate the problem

# simulate data
a <- 3 # intercept
b <- 2 # slope
x <- rnorm(28, 0, 1)
eps <- rnorm(28, 0, 2)
y <- a + b*x + eps  

# put data into list
data_reg <- list(N = 28, x = x, y = y)

# model string with priors sampled directly from generated quantities block
ms <- "
  data {
    int<lower=0> N;
    vector[N] x;
    vector[N] y;
  }
  
  parameters {
    real alpha;
    real beta;
    real<lower=0> sigma;
  }
  
  model {
    vector[N] mu;
    sigma ~ cauchy(0, 2);
    beta ~ normal(0,10);
    alpha ~ normal(0,100);
    for ( i in 1:N ) {
    mu[i] = alpha + beta * x[i];
    }
    y ~ normal(mu, sigma);
  }
  
  generated quantities {
    real <lower=0> sigma_sim = cauchy_rng(0, 2);  
    real beta_sim = normal_rng(0, 10);
    real alpha_sim = normal_rng(0, 100);
  }
"

# now fit the model in stan
fit1 <- stan(model_code = ms,
             data = data_reg,     
             chains = 1,            
             warmup = 1e3,      
             iter = 2e3)

Running this model results in lots of error nessages like

Chain 1: Exception: model6993f7efae5_73687731568b4ebb4a09bd0cbbe77f50_namespace::write_array: sigma_sim is -1.09914, but must be greater than or equal to 0  (in 'model6993f7efae5_73687731568b4ebb4a09bd0cbbe77f50' at line 26)

And ultimately in

error occurred during calling the sampler; sampling not done

However if we simply remove the <lower=0> from the specification of thesigma_sim` variable

ms2 <- "
  data {
    int<lower=0> N;
    vector[N] x;
    vector[N] y;
  }
  
  parameters {
    real alpha;
    real beta;
    real<lower=0> sigma;
  }
  
  model {
    vector[N] mu;
    sigma ~ cauchy(0, 2);
    beta ~ normal(0,10);
    alpha ~ normal(0,100);
    for ( i in 1:N ) {
    mu[i] = alpha + beta * x[i];
    }
    y ~ normal(mu, sigma);
  }
  
  generated quantities {
    real sigma_sim = cauchy_rng(0, 2);  
    real beta_sim = normal_rng(0, 10);
    real alpha_sim = normal_rng(0, 100);
  }
"

And run the model

fit2 <- stan(model_code = ms2,
             data = data_reg,     
             chains = 1,            
             warmup = 1e3,      
             iter = 2e3)

It works, but generates values lower than 0.

print(fit2, pars=c("alpha", "beta", "sigma", "alpha_sim", "beta_sim", "sigma_sim"), prob = c(0.05, 0.95)) 

           mean se_mean    sd      5%    95% n_eff Rhat
alpha      2.40    0.02  0.47    1.62   3.16   961    1
beta       1.58    0.02  0.48    0.84   2.41   829    1
sigma      2.40    0.02  0.37    1.90   3.07   561    1
alpha_sim -1.61    3.14 99.91 -169.58 161.22  1012    1
beta_sim   0.35    0.32 10.10  -16.05  16.78   984    1
sigma_sim -0.83    1.61 51.23  -12.40  14.09  1007    1

So how do I truncate _rng() generators in the generated quantities block?

2

Just apply the fabs function to take the absolute value and “fold” was the samples around zero.

1 Like
3

Thank you @betanalpha, how, or rather where do you do this, inside the model string or outside once the samples have been generated?

4 
generated quantities {
    real <lower=0> sigma_sim = fabs(cauchy_rng(0, 2));   }
2 Likes