# Truncating a cauchy_rng() sampler

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?

Just apply the fabs function to take the absolute value and â€śfoldâ€ť was the samples around zero.

1 Like

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

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