Integrate_1d() error: "The sinh_sinh quadrature evaluated your function at a singular point, leading to the value nan."

[karimn]

Hi,

I previously had a thread on some similar issues with integrate_1d() here.

https://discourse.mc-stan.org/t/calling-integrate-1d-reports-parsing-error-parser-expected/19611

I’m starting a new thread since the previous one was dealing with all kinds of issues as I was learning my way around this functionality. Previously, @bbbales helped me rewrite my code and I was able to fix most of my problems. But I notice that I occasionally get the below message. I’ve been ignoring this for the time being, but I’m now trying to figure out if that messes up my posterior.

Chain 1 Exception: Exception: Exception: Error in function boost::math::quadrature::sinh_sinh<double>::integrate: The sinh_sinh quadrature evaluated your function at a singular point, leading to the value nan. 
Chain 1 Exception: Exception: Exception: Error in function boost::math::quadrature::sinh_sinh<double>::integrate: The sinh_sinh quadrature evaluated your function at a singular point, leading to the value nan.
Chain 1 sinh_sinh quadrature cannot handle singularities in the domain. 
Chain 1 If you are sure your function has no singularities, please submit a bug against boost.math 
Chain 1  (in '/home/karim/Code/takeup/stan_models/takeup_functions.stan', line 88, column 2, included from 
Chain 1 '/tmp/RtmpXSRlSq/model-43e777e2eb735.stan', line 2, column 0) (in '/home/karim/Code/takeup/stan_models/takeup_functions.stan', line 88, column 2, included from 
Chain 1 '/tmp/RtmpXSRlSq/model-43e777e2eb735.stan', line 2, column 0) (in '/home/karim/Code/takeup/stan_models/takeup_functions.stan', line 88, column 2, included from 
Chain 1 '/tmp/RtmpXSRlSq/model-43e777e2eb735.stan', line 2, column 0)

Below is my code (two small functions). I call expected_delta from another function that is called by algebra_solver_newton. expected_delta uses integrate_1d calling expected_delta_part where the problem happens.

real expected_delta_part(real v, real xc, real[] theta, data real[] x_r, data int[] x_i) {
  real w = theta[1];
  real u_sd = theta[2];
  
  real wmv_cdf = Phi_approx((w - v) / u_sd);
  real v_lpdf = normal_lpdf(v | 0, 1);

  return v * exp(v_lpdf) * wmv_cdf; // <---------- This is the line reported as the problem
}

real expected_delta(real w, real total_error_sd, real u_sd, data real[] x_r, data int[] x_i) {
  real delta_part;
  real F_w; 
 
  print("w = ", w, ", u_sd = ", u_sd); 
   
  delta_part = integrate_1d(expected_delta_part, negative_infinity(), positive_infinity(), { w, u_sd }, x_r, x_i, 0.00001);
  
  print("delta_part = ", delta_part);
  
  F_w = Phi_approx(w / total_error_sd);
  
  return - delta_part / (F_w * (1 - F_w));
}

I exported this function using rstan and tested in R’s integrate function and seems to work fine. I tried all kinds of extreme values for the theta variables (w and u_sd) but the function works fine unless I pass it a NaN.

In the above Stan code, I also added a print statement to try and show the values passed into integrate_1d that cause the problem but I can’t seem to catch it. Below is the kind of output I’m seeing (running a single chain). I expected the exception to be raised between printing out w and u_sd and printing delta_part.

Chain 1 w = 0.125534, u_sd = 0.836498 
Chain 1 delta_part = -0.304586 
Chain 1 w = 0.125534, u_sd = 0.836498 
Chain 1 delta_part = -0.304586 
Chain 1 Exception: Exception: Exception: Error in function boost::math::quadrature::sinh_sinh<double>::integrate: The sinh_sinh quadrature evaluated your function at a singular point, leading to the value nan. 
Chain 1 Exception: Exception: Exception: Error in function boost::math::quadrature::sinh_sinh<double>::integrate: The sinh_sinh quadrature evaluated your function at a singular point, leading to the value nan.
Chain 1 sinh_sinh quadrature cannot handle singularities in the domain. 
Chain 1 If you are sure your function has no singularities, please submit a bug against boost.math 
Chain 1  (in '/home/karim/Code/takeup/stan_models/takeup_functions.stan', line 88, column 2, included from 
Chain 1 '/tmp/RtmpXSRlSq/model-43e777e2eb735.stan', line 2, column 0) (in '/home/karim/Code/takeup/stan_models/takeup_functions.stan', line 88, column 2, included from 
Chain 1 '/tmp/RtmpXSRlSq/model-43e777e2eb735.stan', line 2, column 0) (in '/home/karim/Code/takeup/stan_models/takeup_functions.stan', line 88, column 2, included from 
Chain 1 '/tmp/RtmpXSRlSq/model-43e777e2eb735.stan', line 2, column 0) 

Any suggestions on how I can figure out the cause of this exception?

I’m using the R package {cmdstanr} and using cmdstan version 2.26.

Thanks!
Karim

[bbbales2]

So you can run this function at whatever w and u_sd you like with the exported function, but then when it comes time to run the model it blows up?

If this is the case, then it’s probably the gradient calculation that is exploding.

We do Leibniz integral rule - Wikipedia for computing the gradients of the integrals in the ODE. I dunno what the function looks like, but it’s probably possible for a singularity to appear in a derivative that wasn’t there in the function itself. I don’t see where that would happen here tho, so I’m guessing some sorta numerical thing?

The only things I can think of off the top of my head is try to limit the domain of integration. Maybe just do [-1000, 1000]? Or set extreme values in the integrand to zero (like here).

What are valid ranges for w and u_sd by the way?

[karimn]

It only blows up in the model sometimes. I’m trying to figure out when this happens, hence the print statements that don’t seem to work because the exception seems to happen before the values of u_sd and w are printed. Can the exception be thrown somewhere other than the integrate_1d call?

I tried to change the integration to be limited to [-1000, 1000], [-100, 100], [-10, 10], and [-5, 5] but I’m still getting some errors. I also added a check in the function to make sure I don’t try to call Phi_approx with numbers that are too high or too low and setting the value to 0 and 1 manually.

The error did change though to complaining about tanh_sinh instead of sinh_sinh:

Chain 2 Exception: Exception: Exception: Error in function tanh_sinh<double>::integrate: The tanh_sinh quadrature evaluated your function at a singular point and got nan. Please narrow the bounds of integration or check your function for singularities. (in '/home/karim/Code/takeup/stan_models/takeup_functions.stan', line 97, column 2, included from

To give you an idea of what I’m trying to integrate

\begin{equation} \int_{-\infty}^\infty v \cdot \Phi\left(\frac{w - v}{\sigma}\right)\phi(v)\,\mathrm{d}v, \end{equation}

where \sigma is u_sd. Now, the integration is over a restricted range [-5,5] and I’m manually setting the \Phi to 0 and 1 as appropriate if w-v is too big or small.

To answer your question on the values of u_sd and w: u_sd and has a \mathtt{Normal}^+(0,1) prior and w is used in a probit like estimation (with \sigma as the SD) so shouldn’t be too big.

To show you what the integrand function actually looks like here’s a plot, given varying levels of w and u_sd (Some of these values are extreme but I wanted to see what they look like). I’m also displaying the integration result from R’s integrate function.

This is from the below R code after exposing a modified version of the Stan function that I’m trying to integrate (R’s integrate wants a vector to be passed in), also shown below.

expand.grid(
  w = seq(-10, 10, 2),
  u_sd = seq(0.5, 5.0, 1)
) %>% 
  mutate(
    int_y = map2_dbl(w, u_sd, ~ integrate(expected_delta_part, -10, 10, w = ..1, u_sd = ..2)$value),
    y = map2(w, u_sd, ~ tibble(v = ..3, y = expected_delta_part(v, w = ..1, u_sd = ..2)), v = seq(-10, 10, 0.1)),
  ) %>%
  ggplot() +
  geom_line(aes(v, y), data = . %>% unnest(y)) +
  geom_text(aes(0, 0.35, label = sprintf("%.4f", int_y)), color = alpha("red", 0.5), size = 6, fontface = "bold") +
  coord_cartesian(xlim = c(-5, 5)) +
  facet_grid(vars(u_sd), vars(w), labeller = label_both) +
  theme_minimal() +
  theme(
    axis.text.x = element_blank()
  )

vector expected_delta_part(vector v, real w, real u_sd) {
    int n = num_elements(v);
    vector[n] r;
    
    for (i in 1:n) {
      real std_wmv = (w - v[i]) / u_sd; 
      real v_lpdf = normal_lpdf(v[i] | 0, 1);
      real wmv_cdf;
      
      if (std_wmv > 10) { 
        wmv_cdf = 1; 
      } else if (std_wmv < -10) {
        wmv_cdf = 0; 
      } else {
        wmv_cdf = Phi_approx(std_wmv);
      }
      
      r[i] = v[i] * exp(v_lpdf) * wmv_cdf; 
    }
    
    return r;
  }

[karimn]

When I tried to replace the Phi_approx call (in expected_delta_part) with normal_lcdf, I started getting these errors instead:

Chain 2 Exception: Exception: Exception: Exception: normal_lcdf: Random variable is -nan, but must be not nan! (in '/home/karim/Code/takeup/stan_models/takeup_functions.stan', line 82, column 2, included from 

This call is passed in w - v and since v (I presume) is coming from integrate_1d, the problem must be with w being nan, right? But why isn’t Phi_approx complaining about the same problem in the original code?

I’m also occasionally seeing this error. I don’t know what error flag -5 means (the above integration is happening inside an algebra_solver_newton` call).

Chain 4 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 4 Exception: algebra_solver failed with error flag -5. (in '/home/karim/Code/takeup/stan_models/takeup_functions.stan', line 145, column 2, included from
Chain 4 '/home/karim/Code/rtemp/Rtmp6U47cc/model-8d5097fbd00da.stan', line 2, column 0)
Chain 4 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
Chain 4 but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
Chain 4 

[bbbales2]

@karimn

the above integration is happening inside an algebra_solver_newton` call

Does this happen without an algebra_solver_newton call?

This looks pretty complicated. I’ll set up a model later today and see if I can reproduce it in a model without the algebra_solver (if indeed the problem still happens there edit: which it totally might I’m just not sure)

[karimn]

I haven’t tried without algebra_solver. I’ll give that a try on my end too.

[bbbales2]

Here’s a simple test model I was able to run:

functions {
    real expected_delta_part(real v, real xc, real[] theta, data real[] x_r, data int[] x_i) {
    real w = theta[1];
    real u_sd = theta[2];
    
    real wmv_cdf = Phi_approx((w - v) / u_sd);
    real v_lpdf = normal_lpdf(v | 0, 1);

    return v * exp(v_lpdf) * wmv_cdf; // <---------- This is the line reported as the problem
    }

    real expected_delta(real w, real total_error_sd, real u_sd, data real[] x_r, data int[] x_i) {
        real delta_part;
        real F_w; 
        
        print("w = ", w, ", u_sd = ", u_sd); 
        
        delta_part = integrate_1d(expected_delta_part, negative_infinity(), positive_infinity(), { w, u_sd }, x_r, x_i, 0.00001);
        
        print("delta_part = ", delta_part);
        
        F_w = Phi_approx(w / total_error_sd);
        
        return - delta_part / (F_w * (1 - F_w));
    }
}

transformed data {
    real x_r[0];
    int x_i[0];
}

parameters {
  real w;
  real<lower = 0.5> u_sd;
}

transformed parameters {
  real delta = expected_delta(w, 1.0, u_sd, x_r, x_i);
}

model {
  w ~ normal(0, 3.3);
  u_sd ~ normal(0, 1);
}

With those priors (which I picked based on the test range above), there’s lots of big numbers and divergences:

>>> fit.summary()
               Mean          MCSE        StdDev    5%    50%        95%   N_Eff  N_Eff/s  R_hat
name                                                                                           
lp__  -2.100000e+00  2.900000e-02  1.000000e+00 -4.10 -1.800       -1.1  1200.0    280.0    1.0
w     -7.600000e-02  7.700000e-02  3.200000e+00 -5.40  0.045        5.0  1700.0    390.0    1.0
u_sd   1.100000e+00  1.200000e-02  5.100000e-01  0.55  1.000        2.2  2000.0    440.0    1.0
delta  4.200000e+41  4.200000e+41  2.600000e+43  1.00  6.400  2600000.0  4000.0    900.0    1.0

But the model runs. If I make those priors tiny (set both sds in the priors to 0.1), then the divergences go away and the output looks pretty clean:

>>> fit.summary()
         Mean    MCSE  StdDev    5%     50%   95%   N_Eff  N_Eff/s  R_hat
name                                                                     
lp__  -2.2000  0.0300    1.10 -4.30 -1.9000 -1.10  1200.0    370.0    1.0
w     -0.0035  0.0025    0.10 -0.16 -0.0053  0.17  1700.0    500.0    1.0
u_sd   1.1000  0.0110    0.53  0.55  1.0000  2.20  2300.0    690.0    1.0
delta  1.1000  0.0053    0.23  0.67  1.1000  1.40  2000.0    590.0    1.0

So this problem might require the algebra solver and the integerator to get it to go. So I guess what’s the simplest test of that?

Edit: added missing part of model

[karimn]

Ok, I think I found a minimal model to replicate this issue. Notice that in expected_delta_part I’m using normal_lcdf which reports being passed -nan values. If you comment out lines 8 and 19 and uncomment 10-18, you’ll use Phi_approx instead and it will return the sinh_sinh errors.

Straightforward to run:

library(cmdstanr)

test_model <- cmdstan_model("test.stan")

test_model$sample(
  # seed = 1001234565,
  iter_warmup = 1000,
  iter_sampling = 1000, 
  parallel_chains = 4
)

functions {
  real expected_delta_part(real v, real xc, real[] theta, data real[] x_r, data int[] x_i) {
    real w = theta[1];
    real u_sd = theta[2];
   
    real std_wmv = (w - v) / u_sd; 
    real v_lpdf = normal_lpdf(v | 0, 1);
    real wmv_lcdf = normal_lcdf(w - v | 0, u_sd);
    
    // if (std_wmv > 5) {
    //   wmv_cdf = 1;
    // } else if (std_wmv < -5) {
    //   wmv_cdf = 0;
    // } else {
    //   wmv_cdf = Phi_approx(std_wmv);
    // }
    // 
    // return v * exp(v_lpdf) * wmv_cdf;
    return v * exp(v_lpdf + wmv_lcdf);
  }
  
  real expected_delta(real w, real total_error_sd, real u_sd, data real[] x_r, data int[] x_i) {
    real F_w = Phi_approx(w / total_error_sd); 
   
    real delta_part = integrate_1d(expected_delta_part, negative_infinity(), positive_infinity(), { w, u_sd }, x_r, x_i, 0.00001);
    
    return - delta_part / (F_w * (1 - F_w));
  }
  
  vector vec_expected_delta_part(vector v, real w, real u_sd) {
    int n = num_elements(v);
    vector[n] r;

    for (i in 1:n) {
      r[i] = expected_delta_part(v[i], 0.0, { w, u_sd }, { 0.0 }, { 0 });
    }

    return r;
  }
  
  vector w_fixedpoint_solution_normal(vector model_param, vector theta, data real[] x_r, data int[] x_i) {
    real w = model_param[1];
    
    real beta = theta[1];
    real mu = theta[2];
    real u_sd = theta[3];
    
    real delta = expected_delta(w, sqrt(1 + square(u_sd)), u_sd, x_r, x_i);
    
    return [ w + beta + mu * delta ]';
  }
}

transformed data {
  real x_r[1] = { 0.0 };
  int x_i[1] = { 0 };
}

parameters {
  real<lower = 0> u_sd; 
  real w;
  real<lower = 0> mu;
  real beta;
}

transformed parameters {
  real delta;
  real w2;
  
  {
    real total_sd = sqrt(1 + square(u_sd));
    
    delta = expected_delta(w, total_sd, u_sd, x_r, x_i);
    
    w2 = algebra_solver_newton(w_fixedpoint_solution_normal, [ -beta ]', [ beta, mu, u_sd ]', x_r, x_i)[1];
  }
}

model {
  u_sd ~ normal(0, 1);
  w ~ normal(0, 2);
  mu ~ normal(0, 1);
  beta ~ normal(0, 0.1);
}

[karimn]

We must have posted at almost the same time! Yes, calling the integrator by itself works fine. This seems to happen from inside algebra_solver_newton. See above code.

[karimn]

Why is your code so nicely color coded and not mine? :-)

[rok_cesnovar]

For Stan syntax coloring you need to place it in

```stan
```

[bbbales2]

Will probably get around to this tomorrow AM. Thanks for code sample!

[karimn]

Thanks, Ben.

[bbbales2]

I switched from algebra_solver_newton to algebra_solver and got less warnings, but I think the problem is the equation doesn’t always have a solution:

library(rstan)
rstan_options(javascript = FALSE)

expose_stan_functions("model.stan")

x_r = c(1.0, 2.0)
x_i = c(1, 2)

f = Vectorize(function(x) {
  w_fixedpoint_solution_normal(x, c(0.187118, 1.01727, 0.257311), x_r, x_i)
})

curve(f(x), -5, 5)

ALso the printouts were weird. What I did was:

print(beta, " ", mu, " ", u_sd);
w2 = algebra_solver(w_fixedpoint_solution_normal, [ -beta ]', [ beta, mu, u_sd ]', x_r, x_i)[1];
print("asdf");

The idea was the print that caused the problem wouldn’t have a matching “asdf”.

And then the output looked like:

Chain 1 0.0439722 0.847346 0.265992 
Chain 1 asdf 
Chain 1 Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
Chain 1 Exception: ... // Tons more error messages
Chain 1 0.187118 1.01727 0.257311 
Chain 1 0.203986 0.0980935 1.04618 
Chain 1 asdf 

And so I guess the error was printing before the message that should have printed before the error, but the print itself was still there (without the matching “asdf”).

[karimn]

I switched to algebra_solver and I’m not getting the sinh_sinh or -nan to normal_lcdf errors anymore. I also compared R and Stan’s output for the function being used here and everything looks fine. Solutions can be found within the right range of parameters. Here I’m plotting the integration results and solution from w_fixedpoint_solution_normal. They’re identical and there is a solution shown as a vertical dashed line (I’m using a different sent of values for theta than you did and limiting the range of x to -2.5 and 2.5). I can’t do this with algebra_solver_newton because my version of rstan (from CRAN) doesn’t seem to support that function.

rstan::expose_stan_functions("test.stan")

x_r = c(1.0, 2.0)
x_i = c(1, 2)

theta <- c(beta = -0.187118, mu = 0.1727, u_sd = 0.257311)
# theta <- c(beta = 0.187118, mu = 1.01727, u_sd = 0.257311)

tibble(
  x = seq(-2.5, 2.5, 0.1),
  F_w_r = pnorm(x, sd = sqrt(1 + theta["u_sd"]^2)), 
) %>% 
  rowwise() %>%
  mutate(
    int_delta_part_stan = integrate_delta_part(x, theta["u_sd"], x_r, x_i), 
    int_delta_part_r = integrate(vec_expected_delta_part, -Inf, Inf, w = x, u_sd = theta["u_sd"], x_r, x_i)$value,
  
    delta_stan = expected_delta(x, sqrt(theta["u_sd"]^2 + 1), theta["u_sd"], x_r, x_i), 
    delta_r = - int_delta_part_r / (F_w_r * (1 - F_w_r)), 
  
    w_stan = w_fixedpoint_solution_normal(x, theta, x_r, x_i),
  ) %>% 
  ungroup() %>% 
  mutate(
    w_r = x + theta["beta"] + theta["mu"] * delta_r
  ) %>%  
  select(-F_w_r) %>% 
  pivot_longer(-x, names_to = c("name", "src"), names_pattern = "(.+)_(r|stan)") %>% 
  ggplot(aes(x, value)) +
  geom_line(aes(color = src)) +
  geom_vline(aes(xintercept = x), linetype = "dashed", 
             data = tibble(x = find_fixedpoint(theta["beta"], theta["mu"], theta["u_sd"], x_r, x_i), 
                           name = "w")) +
  facet_wrap(vars(name), scales = "free_y") +
  NULL

I also added this function to the model

  real find_fixedpoint(real beta, real mu, real u_sd, data real[] x_r, data int[] x_i) {
    return algebra_solver(w_fixedpoint_solution_normal, [ -beta ]', [ beta, mu, u_sd ]', x_r, x_i)[1];
  }

I had the same problem with print statements before and I thought perhaps the exception was preventing both the before and after output; I wasn’t seeing any -nan values that normal_lcdf was reporting.

Ok, so the solution seems to be to use algebra_solve and not algebra_solver_newton. I’m going to continue to get error messages as the model strays into values of parameters that don’t have a solution. Would this cause problems for my inference?

[bbbales2]

@karimn

Would this cause problems for my inference?

As long as these only happen in warmup (especially very early warmup), I wouldn’t worry too much about it.

If they happen in sampling or later in warmup, then it would probably be worth investigating what is going on.

If this is the case then you may need to use priors on constraints to stay away from these invalid pieces of parameter space.

[karimn]

Great. Thanks for all your help!