Help writing _rng function [error: require unconstrained. found range constraint.]

[ignacio]

I’m trying to write a function to generate data. This is my model:

\begin{aligned} y_{i,s,t} \sim N(a[i] + \gamma y_{i,s,t-1} , \sigma) \\ a[i] \sim N(\mu,sigma_a) \\ \mu \sim N(700, 75) \\ \sigma_a \sim N_+(0, 100) \\ \gamma \sim N_+(0.5,1) \\ \sigma \sim \Gamma(2,8) \end{aligned}

and this is my stan code:

functions {
  matrix parc_rng(int S, vector pre_test) {
    int N = cols(pre_test); 
    matrix[S, N] out; // holds prior predictive draws
    for (s in 1:S) {
      real mu ;
      vector[N] a_std;              
      real<lower=0> sigma_a;
      real<lower=0> gamma;
      real<lower=0> sigma;
      
      vector[N] a_std;
      vector[N] a;
      vector[N] theta;
      
      mu = normal_rng(700,75);
      sigma_a = normal_rng(0,100)
      gamma = normal_rng(0.5,1);
      sigma = gamma_rng(2,8);
      
      for(n in 1:N){
        a_std[n] = normal_rng(0,1);
        a[n] = mu + sigma_a * a_std[n];
        theta[n] = a[n] + gamma*pre_test[n];
        out[s, n] = normal_rng(theta[n], sigma)
      }
      
    return out;
  }
}

This is the error that I get when I try to expose the function to R:

library(rstan)
expose_stan_functions('parcc_rng.stan')

SYNTAX ERROR, MESSAGE(S) FROM PARSER:

require unconstrained. found range constraint.
  error in 'model_parcc_rng' at line 8, column 28
  -------------------------------------------------
     6:       real mu ;
     7:       vector[N] a_std;              
     8:       real<lower=0> sigma_a;
                                   ^
     9:       real<lower=0> gamma;
  -------------------------------------------------

Error in stanc(file = stanmodel, allow_undefined = TRUE, obfuscate_model_name = FALSE) : 
  failed to parse Stan model 'parcc_rng' due to the above error.

If I’m understanding this correctly the problem is that I’m setting a lower bound for some of my parameters. Is that correct? If so, how can I fix my problem?

Thanks!

[shoshievass]

Jim is probably writing a better reply right now, but my guess is tl;dr you can’t define [constrained] params inside a function. All parameter defs (to be treated as params, etc.) only belong in the parameter block w/ priors in the model section.

What you want to do is draw from constrained rng functions, which you might have to define externally, either by clever tricks (see Jim’s post - this will generally work better) or by writing out separate rng functions that throw out draws outside your constraints (there are clever-ish versions and simple versions of this).

[James_Savage]

Just define sigma as unconstrained and use exp(sigma) where you use it. Naturally you’ll not want to do exp(normal_rng(0, 100)).

Or you could guarantee positivity by generating draws from a truncated normal by providing random draws from a uniform

https://gist.github.com/khakieconomics/2c78e26a41474c6c653f35958c0a8e61

[James_Savage]

Alternatively draw prior predictive values from your model with this technique (this has the useful property of testing your model fitting code and will throw warnings if your parameterization is bad)

An easy way to simulate fake data from your Stan model

[ignacio]

@James_Savage

Or you could guarantee positivity by generating draws from a truncated normal by providing random draws from a uniform

Sorry for the dumb question, but could you show me how to use truncnorm_ng. For example, instead of normal_rng(0,100); i should call truncnorm_ng(p, 0, positive_infinity(), 0,100) What should i use for p?

[James_Savage]

you could do something like

for(s in 1:S) {
  real u = uniform_rng(0.0, 1.0);
  real u2 = uniform_rng(0.0, 1.0);
  sigma_a = truncnorm_ng(u, 0, positive_infinity(), 0,100);
  gamma = truncnorm_ng(u2, 0, positive_infinity(), 0.5, 1);
  // ... 
}

[ignacio]

@James_Savage

you could do something like

I’m getting this error:

No matches for: 

  truncnorm_ng(real, int, real, int, int)

Available argument signatures for truncnorm_ng:

  truncnorm_ng(vector, real, real, real, real)

  error in 'model_parcc_rng' at line 57, column 63
  -------------------------------------------------
    55:       real u = uniform_rng(0.0, 1.0);
    56:       real u2 = uniform_rng(0.0, 1.0);
    57:       sigma_a = truncnorm_ng(u, 0, positive_infinity(), 0,100);
                                                                      ^
    58:       gamma = truncnorm_ng(u2, 0, positive_infinity(), 0.5, 1);
  -------------------------------------------------

Error in stanc(file = stanmodel, allow_undefined = TRUE, obfuscate_model_name = FALSE) : 
  failed to parse Stan model 'parcc_rng' due to the above error.

https://gist.github.com/ignacio82/855ac1d87d34e4ffc5d5bf0240e63b59

[James_Savage]

Ug sorry – that rng is vectorized in p. You’ll have to change the signature of p to be real, alter the for loop, and don’t use integers when you call it (ie use 0.0, 1.0 rather than 0, 1 respectively)

[ignacio]

Thanks @James_Savage ! I’m almost there. The problem is that my code is producing a lot of NaN for s>1 This is my stan code:

functions {
  
  // lower bound is a, upper bound is b, rv is x, mean is mu, sd is sigma
  vector xi(vector x, real mu, real sigma) {
    return((x - mu)./sigma);
  }
  real alpha(real a, real mu, real sigma) {
    real out;
    out = (a==negative_infinity())? negative_infinity(): (a - mu)/sigma;
    return(out);
  }
  real beta(real b, real mu, real sigma) {
    real out;
    out = (b==positive_infinity())? positive_infinity(): (b - mu)/sigma;
    return(out);
  }
  real Z(real a, real b, real mu, real sigma) {
    return(normal_cdf(beta(b, mu, sigma), 0.0, 1.0) - normal_cdf(alpha(a, mu, sigma), 0.0, 1.0));
  }
  
  // converts a vector of inverse quantiles (a CDF) to the corresponding quantile of 
  // a truncated normal distribution in [a, b] with location = location and scale = scale. 
  // convenient for 
  real truncnorm_ng(real p, // vector the cdf of some random variable
                      real a, // lower bound, can be negative_infinity() and positive_infinity()
                      real b, // upper bound, can be negative_infinity() and positive_infinity()
                      real location, // location and scale are self explanatory
                      real scale) {
    real out;
    real tmp_Z;
    real tmp_alpha;
    
    tmp_alpha = normal_cdf(alpha(a, location, scale), 0, 1);
    tmp_Z = normal_cdf(beta(b, location, scale), 0, 1) - tmp_alpha;
    out = inv_Phi(tmp_alpha + p*tmp_Z)*scale + location;
    return(out);
  }
  
  matrix parcc_rng(int S, vector pre_test) {
    int N = rows(pre_test); 
    matrix[S, N] out; // holds prior predictive draws
    for (s in 1:S) {
      real mu ;
      real sigma_a;
      real gamma;
      real sigma;
      
      vector[N] a_std;
      vector[N] a;
      vector[N] theta;
      
      real u = uniform_rng(0.0, 1.0);
      real u2 = uniform_rng(0.0, 1.0);
      sigma_a = truncnorm_ng(u, 0, positive_infinity(), 0,100);
      gamma = truncnorm_ng(u2, 0, positive_infinity(), 0.5, 1);
      
      mu = normal_rng(700,75);
      sigma = gamma_rng(2,8);
      
      for(n in 1:N){
        a_std[n] = normal_rng(0,1);
        a[n] = mu + sigma_a * a_std[n];
        theta[n] = a[n] + gamma*pre_test[n];
        out[s, n] = normal_rng(theta[n], sigma);
      }
      
    return out;
  }
}
}

and this is my r code:

library(rstan)
expose_stan_functions('parcc_rng.stan')

test <- parcc_rng(S = 10, pre_test = c(700,750), seed = 9782L)
test

          [,1]     [,2]
 [1,] 1496.387 1778.769
 [2,]      NaN      NaN
 [3,]      NaN      NaN
 [4,]      NaN      NaN
 [5,]      NaN      NaN
 [6,]      NaN      NaN
 [7,]      NaN      NaN
 [8,]      NaN      NaN
 [9,]      NaN      NaN
[10,]      NaN      NaN

What am I missing? :-(

[James_Savage]

Look at the 4th last line. Pretty sure you don’t want that inside your for loop.

The thing I’m concerned about is why you’re not saving the parameter draws for each student. The whole point of simulating data is to see how your model goes at recapturing the “known unknowns”. I strongly recommend you use the second method I recommended in my blog, as it’ll save all the parameters, transformed parameters, and simulated data for each draw.

[ignacio]

@James_Savage

Look at the 4th last line. Pretty sure you don’t want that inside your for loop.

Thanks! That was a not very smart mistake on my part :(

@James_Savage

The thing I’m concerned about is why you’re not saving the parameter draws for each student. The whole point of simulating data is to see how your model goes at recapturing the “known unknowns”. I strongly recommend you use the second method I recommended in my blog, as it’ll save all the parameters, transformed parameters, and simulated data for each draw.

This is the next thing on my to-do list.

I’m leaving a working version of my code here in case is useful for someone else in the future:

functions {
  
  // lower bound is a, upper bound is b, rv is x, mean is mu, sd is sigma
  vector xi(vector x, real mu, real sigma) {
    return((x - mu)./sigma);
  }
  real alpha(real a, real mu, real sigma) {
    real out;
    out = (a==negative_infinity())? negative_infinity(): (a - mu)/sigma;
    return(out);
  }
  real beta(real b, real mu, real sigma) {
    real out;
    out = (b==positive_infinity())? positive_infinity(): (b - mu)/sigma;
    return(out);
  }
  real Z(real a, real b, real mu, real sigma) {
    return(normal_cdf(beta(b, mu, sigma), 0.0, 1.0) - normal_cdf(alpha(a, mu, sigma), 0.0, 1.0));
  }
  
  // converts a vector of inverse quantiles (a CDF) to the corresponding quantile of 
  // a truncated normal distribution in [a, b] with location = location and scale = scale. 
  // convenient for 
  real truncnorm_ng(real p, // vector the cdf of some random variable
                      real a, // lower bound, can be negative_infinity() and positive_infinity()
                      real b, // upper bound, can be negative_infinity() and positive_infinity()
                      real location, // location and scale are self explanatory
                      real scale) {
    real out;
    real tmp_Z;
    real tmp_alpha;
    
    tmp_alpha = normal_cdf(alpha(a, location, scale), 0, 1);
    tmp_Z = normal_cdf(beta(b, location, scale), 0, 1) - tmp_alpha;
    out = inv_Phi(tmp_alpha + p*tmp_Z)*scale + location;
    return(out);
  }
  
  matrix parcc_rng(int S, vector pre_test) {
    int N = rows(pre_test); 
    matrix[S, N] out; // holds prior predictive draws
    for (s in 1:S) {
      real mu ;
      real sigma_a;
      real gamma;
      real sigma;
      
      vector[N] a_std;
      vector[N] a;
      vector[N] theta;
      
      real u = uniform_rng(0.0, 1.0);
      real u2 = uniform_rng(0.0, 1.0);
      sigma_a = truncnorm_ng(u, 0, positive_infinity(), 0,100);
      gamma = truncnorm_ng(u2, 0, positive_infinity(), 0.5, 1);
      
      mu = normal_rng(700,75);
      sigma = gamma_rng(2,8);
      
      for(n in 1:N){
        a_std[n] = normal_rng(0,1);
        a[n] = mu + sigma_a * a_std[n];
        theta[n] = a[n] + gamma*pre_test[n];
        out[s, n] = normal_rng(theta[n], sigma);
      }
  }
  return out;
}
}