# Problem on how to use Integrate_1d function

I want to use integrate_1d function to calculate the integral of a complex function. But I met an error when parsing it.

Here is the error:
Ill-typed arguments supplied to function ‘integrate_1d’. Available signatures:
((real, real, real[], data real[], data int[]) => real, real, real, real[], data real[], data int[]) => real
((real, real, real[], data real[], data int[]) => real, real, real, real[], data real[], data int[], data real) => real
Instead supplied arguments of incompatible type: (real, real, real[], real[], int[]) => real, real, real, real, real[,], int[].

This is my code. I think I have already put x_i into transformed data block, and use the variable in data for x_r. But I have no idea why there is such error. Is there any idea how to solve this problem?

functions {
real normal_density_nc(real x,          // Function argument
real xc,         // Complement of function argument
//  on the domain (defined later)
real[] theta,    // parameters
real[] x_r,      // data (real)
int[] x_i) {     // data (integer)
real mu = x_r[1];
real sigma = x_r[2];

return 1 / (sqrt(2 * pi()) * sigma) * exp(-0.5 * ((x - mu) / sigma)^2) * (-1/sigma + (1 / sigma ^ 2) * (x - mu)^2 );

}
}
data {

// Gaussian model
int W;                                                          // Number of weeks (typically 12)
int N;                                                          // Number of subjects
int Xdim;                                                       // Dimension of X - latent low dimensional structure of the phenotype

// Exogenous variables
int exo_q_num;                                                  // number of exogenous survey questions
real U[N,W,exo_q_num];                                       // exogenous survey questions - missing weeks were linearly interpolated outside of Stan

// Change detection
int<lower=1> W_nc_obs[N];                                          // Number of weeks WITH data
int<lower=0> W_nc_mis[N];                                          // Number of weeks WITHOUT data
int<lower=0> idx_nc_obs[N,W];                              // Indices of weeks WITH data
int<lower=1> P_nc;                                              // Number of parameters
int<lower=0> T_max_nc;                                          // Max number of trials across subjects across weeks
int<lower=0> Tr_nc[N, W];                                // Number of trials for each subj for each week
int deltaM[N, W, T_max_nc];                              // Mu of dots - subj x weeks x trials
real TotalS[N, W, T_max_nc];                             // Sd of dots - subj x weeks x trials
int choice_nc[N, W, T_max_nc];                           // choice nc - subj x weeks x trials

real reward[N, W];
int x_i[0];
}

//transformed data{

//}

parameters {
matrix[Xdim, exo_q_num] B;
real weber_nc_pr_std[N,W];                  // nc
matrix[1,Xdim] C;
matrix[N, Xdim] X1;
real<lower=0, upper=1> eta[N]; //eta
vector[1] mu_d;                     // constant offset
vector<lower=0>[1] sigma_r;
}

model {
to_vector(B) ~ normal(0, 0.1);
to_vector(C) ~ normal(0, 0.1);
to_vector(X1) ~ normal(0, 0.1);
mu_d ~ normal(0, 0.1);
sigma_r ~ normal(0, 0.1);
eta ~ normal(0.01, 0.1);// prior for eta

real weber_nc_pr[N,W];           // cd
real weber_nc[N,W];

real X[N,W,Xdim];

real reward_sum;

for (s in 1:N) {
reward_sum = 0;
int w1 = 1;
real weeks =0;
for (w in 1:W) {

reward_sum += reward[s, w];

if (idx_nc_obs[s,w] != 0) {
weeks = weeks + 1;
}

weber_nc_pr_std[s,w] ~ std_normal();

if (w == 1) {
X[s,w,] = to_array_1d(inv_logit(X1[s,]));
} else {
X[s,w,] = to_array_1d(inv_logit(to_vector(X[s,w1,]) + eta[s].*(reward[s, w] - reward_sum/(weeks + 0.0000001))* to_matrix(C)' * to_vector(nc_grad_sum[s, w1])));
}

weber_nc_pr[s,w] =  to_array_1d(mu_d[1] + B * to_vector(U[s,w,]) + C[1,] * to_vector(X[s,w,]) + sigma_r[1] * weber_nc_pr_std[s, w])[1];
weber_nc[s,w]= exp(weber_nc_pr[s,w]);          // nc

if (idx_nc_obs[s,w] != 0) {                                  // if week exists

vector[Tr_nc[s,w]] z = - to_vector(deltaM[s,w,:Tr_nc[s,w]]) ./ (weber_nc[s,w] * to_vector(TotalS[s,w,:Tr_nc[s,w]]));
choice_nc[s,w,:Tr_nc[s,w]] ~ bernoulli(Phi(z));

vector[Tr_nc[s,w]] itc_p = Phi(z);
vector[Tr_nc[s,w]] log_likelihood_p_grad_nc = to_vector(choice_nc[s,w,:Tr_nc[s,w]]) ./ (itc_p + 0.00001) - (1 - to_vector(choice_nc[s,w,:Tr_nc[s,w]])) ./ (1 - itc_p + 0.00001);

real theta[0];

for(t in 1: Tr_nc[s,w]){
negative_infinity(),
0.0,
weber_nc[s, w],
{ deltaM[s, w, t],  TotalS[s, w, t] },
x_i);
}

w1 = w;

}

}
}
}

1 Like

Hi and sorry about the delay. This is not my area but I searched around on the forums here and found something similar here Parametrised integrals

At the very least this will get bumped up so folks who work with integrate_1d might see it!

1 Like

I believe the issue is that the integrand function (normal_density_nc) is expecting the theta input to be of type real[] (an array of reals), but you’re passing a single real (weber_nc[s, w]).

Try adding curly braces around weber_nc[s, w] in your integrate_1d call:

nc_grad[t] = integrate_1d(normal_density_nc,
negative_infinity(),
0.0,
{ weber_nc[s, w] },
{ deltaM[s, w, t],  TotalS[s, w, t] },
x_i);

2 Likes

Thank you so much! But after fixing it, I met several new problems. I totally have no idea about what such error represents. Hope someone could help me!

/Users/eva/Desktop/computational_phenotyping/nc_gradient_policy.hpp:1538:14: note: in instantiation of function template specialization 'nc_gradient_policy_model_namespace::nc_gradient_policy_model::log_prob_impl<true, true, std::__1::vector<stan::math::var_value<double, void>, std::__1::allocator<stan::math::var_value<double, void> > >, std::1::vector<int, std::1::allocator >, nullptr, nullptr>’ requested here
return log_prob_impl<propto
, jacobian
>(params_r, params_i, pstream);
^
/Users/eva/.cmdstan/cmdstan-2.26.0/stan/src/stan/model/model_base_crtp.hpp:191:50: note: in instantiation of function template specialization ‘nc_gradient_policy_model_namespace::nc_gradient_policy_model::log_prob<true, true, stan::math::var_value<double, void> >’ requested here
return static_cast<const M*>(this)->template log_prob<true, true>(
^
^
/Users/eva/Desktop/computational_phenotyping/nc_gradient_policy.hpp:984:24: note: insert an explicit cast to silence this issue
deltaM[(s - 1)][(w - 1)][(t - 1)],
^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
static_cast( )

That doesn’t look to contain the actual error, are you able to post the complete log? If it’s easier, you can copy it into a text file and upload that

1 Like

compile log.txt (29.7 KB)
Here is th compile log. Thank you so much for your help!

I can see the error but I’m not sure where it’s coming from. Would you be able to upload your Stan model?

Sure! Here is the stan model! Hope it helps, thank you so much!

Thanks for that, unfortunately that model doesn’t parse correctly because you have mismatched curly brackets in the model block. Can you fix those up, and post the model again if it still has the compilation error?

Hi, sorry for the late reply. Here is the correct version of the model and its corresponding error. Could you help me to see it? Thank you so much for your help!

The error here is caused by this section:

nc_grad[1, t] = integrate_1d(normal_density_nc,
negative_infinity(),
0.0,
{weber_nc_pr_std[s, w]},
{ deltaM[s, w, t], TotalS[s,w,t] }, x_i
);


Specifically the array: { deltaM[s, w, t], TotalS[s,w,t] }

Your integrating function normal_density_nc expects this to be an array of reals (real[] x_r), but deltaM is an int while TotalS is a real. This causes the c++ to fail (in a very uninformative way).

You need to either declare deltaM as an array of real:

    real deltaM[N, W, T_max_nc];


Or pass it in using the x_i argument.

It works! Thank you so much!! Another question is that is the integrate_1d function really slow to calculate the results? I have run the model but it seems it is really slow…Is there any method to increase the speed of calculating the integral? Thanks!

There are some optimisations for your integrating function that are possible.

1 / (sqrt(2 * pi()) * sigma) * exp(-0.5 * ((x - mu) / sigma)^2) * (-1/sigma + (1 / sigma ^ 2) * (x - mu)^2 ) * 2 * theta[1] * x_r[2] ^ 2;


Which is:

\frac{1}{\sqrt{2\pi}\sigma}\exp\left( -\frac{1}{2}\left( \frac{x - \mu}{\sigma} \right)^2 \right) \cdot \left( - \frac{1}{\sigma} + \frac{1}{\sigma^2} \cdot (x-\mu)^2 \right) \cdot 2 \cdot \theta_1 \cdot x_2^2

The first half of that function is just the PDF of the normal distribution, so we can simplify that to:

\text{Normal}(x|\mu,\sigma) \cdot \left( - \frac{1}{\sigma} + \frac{1}{\sigma^2} \cdot (x-\mu)^2 \right) \cdot 2 \cdot \theta_1 \cdot x_2^2

Then we can move to the log-scale (enabling direct use of the normal_lpdf function):

\log(\text{Normal}(x|\mu,\sigma)) + \log \left( - \frac{1}{\sigma} + \frac{1}{\sigma^2} \cdot (x-\mu)^2 \right) + \log 2 + \log \theta_1 + 2\log \text{xr_2}

And then exponentiate the result before returning

Given that \frac{1}{\sigma} is used twice, we can compute that once and then just re-use it. log(2) is provided by Stan as the constant log2(), so that doesn’t require computing. Finally, given that xr_2 is data, we can pre-compute 2\log \text{xr_2} in the transformed data block and then pass it to the integrating function. The aim is to avoid as much redundant computation as possible.

Putting that all together, your model looks like:

functions {
real normal_density_nc(real x,          // Function argument
real xc,         // Complement of function argument
//  on the domain (defined later)
real[] theta,    // parameters
real[] x_r,      // data (real)
int[] x_i) {     // data (integer)
real mu = x_r[1];
real sigma = x_r[2]* theta[1];
real inv_sigma = inv(sigma);
real x_r2_logsquared = x_r[3];

return exp(normal_lpdf(x | mu, sigma)
+ log(inv_sigma - inv_sigma^2 * (x - mu)^2)
+ log2()
+ log(theta[1])
+ x_r2_logsquared);

}
}
data {

// Gaussian model
int W;                                                          // Number of weeks (typically 12)
int N;                                                          // Number of subjects
int Xdim;                                                       // Dimension of X - latent low dimensional structure of the phenotype

// Exogenous variables
int exo_q_num;                                                  // number of exogenous survey questions
real U[N,W,exo_q_num];                                       // exogenous survey questions - missing weeks were linearly interpolated outside of Stan

// Change detection
int<lower=1> W_nc_obs[N];                                          // Number of weeks WITH data
int<lower=0> W_nc_mis[N];                                          // Number of weeks WITHOUT data
int<lower=0> idx_nc_obs[N,W];                              // Indices of weeks WITH data
int<lower=1> P_nc;                                              // Number of parameters
int<lower=0> T_max_nc;                                          // Max number of trials across subjects across weeks
int<lower=0> Tr_nc[N, W];                                // Number of trials for each subj for each week
int deltaM[N, W, T_max_nc];                              // Mu of dots - subj x weeks x trials
matrix[W, T_max_nc] TotalS[N];                             // Sd of dots - subj x weeks x trials
int choice_nc[N, W, T_max_nc];                           // choice nc - subj x weeks x trials

real reward[N, W];
int x_i[0];
}

transformed data{
matrix[W, T_max_nc] TotalS_logsquared[N];
for(n in 1:N) {
TotalS_logsquared[n] = 2 * log(TotalS[n]);
}
}

parameters {
matrix[Xdim, exo_q_num] B;
real weber_nc_pr_std[N,W];                  // nc
matrix[1,Xdim] C;
matrix[N, Xdim] X1;
real<lower=0, upper=1> eta[N]; //eta
vector[1] mu_d;                     // constant offset
vector<lower=0>[1] sigma_r;
}

model {
to_vector(B) ~ normal(0, 0.1);
to_vector(C) ~ normal(0, 0.1);
to_vector(X1) ~ normal(0, 0.1);
mu_d ~ normal(0, 0.1);
sigma_r ~ normal(0, 0.1);
eta ~ normal(0.01, 0.1);// prior for eta

real weber_nc_pr[N,W];           // cd
real weber_nc[N,W];

real X[N,W,Xdim];

real reward_sum;

for (s in 1:N) {
reward_sum = 0;
int w1 = 1;
real weeks =0;
for (w in 1:W) {

reward_sum += reward[s, w];

if (idx_nc_obs[s,w] != 0) {
weeks = weeks + 1;
}

weber_nc_pr_std[s,w] ~ std_normal();

if (w == 1) {
X[s,w,] = to_array_1d(inv_logit(X1[s,]));
} else {
X[s,w,] = to_array_1d(inv_logit(to_vector(X[s,w1,]) + eta[s].*(reward[s, w] - reward_sum/(weeks + 0.0000001))* to_matrix(C)' * to_vector(nc_grad_sum[s, w1])));
}

weber_nc_pr[s,w] =  to_array_1d(mu_d[1] + B * to_vector(U[s,w,]) + C[1,] * to_vector(X[s,w,]) + sigma_r[1] * weber_nc_pr_std[s, w])[1];
weber_nc[s,w]= exp(weber_nc_pr[s,w]);          // nc

if (idx_nc_obs[s,w] != 0) {                                  // if week exists

vector[Tr_nc[s,w]] z = - to_vector(deltaM[s,w,:Tr_nc[s,w]]) ./ (weber_nc[s,w] * to_vector(TotalS[s,w,:Tr_nc[s,w]]));
choice_nc[s,w,:Tr_nc[s,w]] ~ bernoulli(Phi(z));

vector[Tr_nc[s,w]] itc_p = Phi(z);
vector[Tr_nc[s,w]] log_likelihood_p_grad_nc = to_vector(choice_nc[s,w,:Tr_nc[s,w]]) ./ (itc_p + 0.00001) - (1 - to_vector(choice_nc[s,w,:Tr_nc[s,w]])) ./ (1 - itc_p + 0.00001);

//negative_infinity(),
//0,
// { deltaM[s,w,:Tr_nc[s,w]], weber_nc[s,w] * to_vector(TotalS[s,w,:Tr_nc[s,w]]) }, x_r, x_i);

real theta[0];

for(t in 1: Tr_nc[s,w]){
negative_infinity(),
0.0,
{weber_nc_pr_std[s, w]},
{ deltaM[s, w, t], TotalS[s,w,t], TotalS_logsquared[s,w,t]}, x_i
);
}

w1 = w;

}

}
}
}

2 Likes

Hi, thank you so much for such detailed explanation!

But there was a tiny bug in my function and it should be

functions {
real normal_density_nc(real x,          // Function argument
real xc,         // Complement of function argument
//  on the domain (defined later)
real[] theta,    // parameters
real[] x_r,      // data (real)
int[] x_i) {     // data (integer)
real mu = x_i[1];
real sigma = x_r[1]* theta[1];

return 1 / (sqrt(2 * pi()) * sigma) * exp(-0.5 * ((x - mu) / sigma)^2) * (-1/sigma + (1 / sigma ^ 3) * (x - mu)^2 );

}
}


And then I tried using your method to simplify it by

  return exp(normal_lpdf(x | mu, sigma)
+ log(-1 + inv_sigma ^ 2 * (x - mu)^2)
-log(sigma)
);


But there are very serious overflow problem.

Exception: Error in function boost::math::quadrature::exp_sinh<double>::integrate: The exp_sinh quadrature evaluated your function at a singular point and returned nan. Please ensure your function evaluates to a finite number over its entire domain. (in '/Users/chenwenqi/Desktop/computational_phenotyping/nc_gradient_policy.stan', line 17, column 2 to line 20, column 15)


And the log(-1 + inv_sigma ^ 2 * (x - mu)^2) will be nan at most time. I have also tried the -log(sigma ^2) or directly divide sigma inside but all of them did not work. Is there any better method to solve it? Thank you so much!!

Another thing I need your help is whether such function can be used when running in chains. Actually it returns error when I transformed it into multiprocessing version. I guess it is because the integrate_1d function requires the data parameters are directly from data block but not from another parameter in the function. Is there any method to solve it?

My code is here

functions {
real normal_density_nc(real x,          // Function argument
real xc,         // Complement of function argument
//  on the domain (defined later)
real[] theta,    // parameters
real[] x_r,      // data (real)
int[] x_i) {     // data (integer)
real mu = x_i[1];
real sigma = x_r[1]* theta[1];
return exp(normal_lpdf(x | mu, sigma))  * (-1/sigma + (1 / sigma ^ 3) * (x - mu)^2 );

}

real partial_sum(real[,,] all_parameters_sliced, int start, int end,  matrix B, matrix C, real[,,] U,  matrix X1,
vector sigma_r, vector mu_d, int W, int Xdim, real[] eta,
int[,] idx_nc_obs, int[,,] choice_ncl, int[,,] deltaMl, real[,,] TotalSl,  int[,] Tr_ncl, real[,] reward
){
real lt = 0;

for (n in 1:(end-start+1)) {
real weber_nc_prl[W] = all_parameters_sliced[n,,1];
real weber_ncl[W] = all_parameters_sliced[n,,2];
real X[W,Xdim];

int s = start + (n - 1);

real reward_sum = 0;
int w1 = 1;
real weeks =0;

for (w in 1:W) {
reward_sum += reward[s,w];
if (idx_nc_obs[s,w] != 0) {
weeks = weeks + 1;
}

if (w == 1) {
X[w,] = to_array_1d(inv_logit(X1[s,]));
} else {
X[w,] = to_array_1d(inv_logit(to_vector(X[w1,]) + eta[s] * (reward[s, w] - reward_sum/(weeks+ 0.00000001)) * C' * to_vector(nc_grad_sum[w1]))); // 68
}

lt += normal_lpdf(weber_nc_prl[w] | mu_d[1] + B * to_vector(U[s,w,]) + C[1,] * to_vector(X[w,]),sigma_r[1]);
if (idx_nc_obs[s,w] != 0) {                                  // if week exists

vector[Tr_ncl[s,w]] z = - to_vector(deltaMl[s,w,:Tr_ncl[s,w]]) ./ (weber_ncl[w] * to_vector(TotalSl[s,w,:Tr_ncl[s,w]]));
lt += bernoulli_lpmf(choice_ncl[s,w,:Tr_ncl[s,w]] | Phi(z));

vector[Tr_ncl[s,w]] itc_p = Phi(z);
vector[Tr_ncl[s,w]] log_likelihood_p_grad_nc = to_vector(choice_ncl[s,w,:Tr_ncl[s,w]]) ./ (itc_p + 0.00001) - (1 - to_vector(choice_ncl[s,w,:Tr_ncl[s,w]])) ./ (1 - itc_p + 0.00001);

for(t in 1: Tr_ncl[s,w]){
real nc_integral = integrate_1d(normal_density_nc,
negative_infinity(),
0.0,
{weber_ncl[ w]},
{ TotalSl[s,w,t] }, {deltaMl[s, w, t]}
);

}

w1 = w;

}// if for idx_nc_obs
}// end loop for w

}//end loop for n
}//end partial sum
}//end functions block
data {

// Gaussian model
int W;                                                          // Number of weeks (typically 12)
int N;                                                          // Number of subjects
int Xdim;                                                       // Dimension of X - latent low dimensional structure of the phenotype

// Exogenous variables
int exo_q_num;                                                  // number of exogenous survey questions
real U[N,W,exo_q_num];                                       // exogenous survey questions - missing weeks were linearly interpolated outside of Stan

// Change detection
int<lower=1> W_nc_obs[N];                                          // Number of weeks WITH data
int<lower=0> W_nc_mis[N];                                          // Number of weeks WITHOUT data
int<lower=0> idx_nc_obs[N,W];                              // Indices of weeks WITH data
int<lower=1> P_nc;                                              // Number of parameters
int<lower=0> T_max_nc;                                          // Max number of trials across subjects across weeks
int<lower=0> Tr_nc[N, W];                                // Number of trials for each subj for each week
int deltaM[N, W, T_max_nc];                              // Mu of dots - subj x weeks x trials
real TotalS[N, W, T_max_nc];                             // Sd of dots - subj x weeks x trials
int choice_nc[N, W, T_max_nc];                           // choice nc - subj x weeks x trials

real reward[N, W];
//int x_i[0];
}

parameters {
matrix[Xdim, exo_q_num] B;
real weber_nc_pr[N,W];                  // nc

matrix[1,Xdim] C;
matrix[N, Xdim] X1;
real<lower=0, upper=1> eta[N]; //eta
vector[1] mu_d;                     // constant offset
vector<lower=0>[1] sigma_r;
}

transformed parameters {
real<lower=0> weber_nc[N,W] = exp(weber_nc_pr);
}

model {
to_vector(B) ~ normal(0, 0.1);
to_vector(C) ~ normal(0, 0.1);
to_vector(X1) ~ normal(0, 0.1);
mu_d ~ normal(0, 0.1);
sigma_r ~ normal(0, 0.1);
eta ~ normal(0.01, 0.1);// prior for eta

//real weber_nc_pr[N,W];           // nc
//real weber_nc[N,W];
real all_parameters[N, W, 12];

for (n in 1:N) {
for (w in 1:W) {
all_parameters[n, w, 1] = weber_nc_pr[n,w];
all_parameters[n, w, 2] = weber_nc[n,w];

}
}

target += reduce_sum(partial_sum, all_parameters, 1, B, C, U, X1, sigma_r, mu_d, W, Xdim, eta,
idx_nc_obs, choice_nc, deltaM, TotalS,  Tr_nc, reward);

}


You can mark the function parameter with the data keyword.

real partial_sum(real[,,] all_parameters_sliced, int start, int end,
matrix B, matrix C, real[,,] U,  matrix X1,
vector sigma_r, vector mu_d, int W, int Xdim, real[] eta,
int[,] idx_nc_obs, int[,,] choice_ncl, int[,,] deltaMl,
data real[,,] TotalSl, int[,] Tr_ncl, real[,] reward
){


But it seems Stan can not identify such keyword in the function. After I added “data”, the error returns as “Function bodies must contain a return statement of correct type in every branch.”

That has nothing to do with the keyword, it’s just that your partial_sum never returns anything. I assume you meant to have return lt; at the end of the function.