Normalising constant integral in Stan

Modeling
1

I would like to compute the following power prior in Stan, but I am unsure if I am calculating the integral correctly or if it is even possible in Stan. The integral below is obtainable in closed form, but I would like to calculate it numerically in Stan as I would like the code to work for integrals that cannot be written in closed form.

The power prior, in my example, can be written as

\pi(\lambda, \alpha_0 \mid D_0) \propto \dfrac{\big[\binom{N_1}{n_1}\lambda^{n_1}(1-\lambda)^{N_1 - n_1}\big]^{\alpha_0}\lambda^{a-1}(1-\lambda)^{b-1}}{\int\big[\binom{N_1}{n_1}\lambda^{n_1}(1-\lambda)^{N_1 - n_1}\big]^{\alpha_0}\lambda^{a-1}(1-\lambda)^{b-1} d\lambda}\pi(\alpha_0),

where D_0 = \{N_1, n_1\} is data, \lambda and \alpha_0 are parameters, a and b are data used to construct the prior for \lambda, and \pi(\alpha_0) is a prior on the power parameter.

The integral can be written as

\binom{N_1}{n_1}^{\alpha_0}\int\lambda^{\alpha_0n_1+a-1}(1-\lambda)^{\alpha_0(N_1 - n_1)+b-1}.

The Stan code is shown below:

functions {
  real Integrand(real lambda, real xc, real[] theta, real[] x_r, int[] x_i){
    
    //real lambda = theta[1];
    real beta = theta[1];
    //real a_lambda= x_r[1];
    //real b_lambda= x_r[2];
    //int N = x_i[1];
    //int n_min = x_i[2];
    //int n = x_i[3];
    
    real z = lambda^(n*alpha+a_lambda-1)*(1-lambda)^(alpha*(N-n)+b_lambda-1);
    
    return(z);
  }
}
data {
  int<lower=0> N;
  int<lower=0> n_min;
  int<lower=n_min-1> n;
  real<lower=0> a_lambda;
  real<lower=0> b_lambda;
  real<lower=0> a_alpha;
  real<lower=0> b_alpha;
}
parameters{
  real<lower=0,upper=1> lambda;
  real<lower=0,upper=1> alpha;
}
model{
  real I;
  I = choose(N, n)^alpha*integrate_1d(Integrand, 0, 1, alpha, {a_lambda, b_lambda}, {N,n_min,n}, 1e-8);
  target += alpha*binomial_lpmf(n| N, lambda) - log(I);
  //priors
  lambda ~ beta(a_lambda,b_lambda);
  alpha ~ beta(a_alpha,b_alpha);
}

I am unsure if this is correct as the first argument in the Integrand function, which is the function argument, is a parameter in my example, so should I pass both {lambda, alpha} as parameters in the integrate_1d function? Or just alpha as lambda is the function argument?

I get an error when running the model Identifier 'n' not in scope.

What does this error mean and is there a way to perform the integral in the above model?

2

I got the model to run using the following script

functions {
  real Integrand(real lambda, real xc, real[] theta, real[] x_r, int[] x_i){
    
    real alpha = theta[1];
    real a_lambda= x_r[1];
    real b_lambda= x_r[2];
    int N = x_i[1];
    int n_min = x_i[2];
    int n = x_i[3];
    
    real z = lambda^(n*alpha+a_lambda-1)*(1-lambda)^(alpha*(N-n)+b_lambda-1);
    
    return(z);
  }
}
data {
  int<lower=0> N;
  int<lower=0> n_min;
  int<lower=n_min-1> n;
  real<lower=0> a_lambda;
  real<lower=0> b_lambda;
  real<lower=0> a_alpha;
  real<lower=0> b_alpha;
}
parameters{
  real<lower=0,upper=1> lambda;
  real<lower=0,upper=1> alpha;
}
model{
  real I;
  I = choose(N, n)^alpha*integrate_1d(Integrand, 0, 1, {alpha}, {a_lambda, b_lambda}, {N,n_min,n}, 1e-8);
  target += alpha*binomial_lpmf(n| N, lambda) - log(I);
  //priors
  lambda ~ beta(a_lambda,b_lambda);
  alpha ~ beta(a_alpha,b_alpha);
}

I hope it is doing what I think it is.

3

Just to answer the first question, it means that variable n used in the definition of Integrand is not defined I the scope of the body of the Integrand function. The functions can’t see the data definitions, so you are going to have to pass n as an argument.