# Numerical issue with computing the integral of an exponential spline

Hello,

I want to fit a (IRT-type) model involving a custom probability function. For this custom probability function
(name: espl_lpdf), I need to compute the integral of a function of type exp(c*bspline(x)), wherein bspline(x) is a B-spline with correspondingly defined knots.

When using the model below on test_data, I get error messages refering to a too large error estimate of the integral. I have difficulties in coming up with a proper solution to this problem. After all, the integral I compute (via the integrate_1d function - using the quadrature rule) is well-defined and quadrature should not have any problems with this integral (the function is smooth - except for the knot-locations; the integral is one-dimensional).

I would be very grateful for any comments, advices on this topic and for suggestions, on how to resolve this problem. For the sake of clarity, I have appended the whole Stan- and R-code – as well as the corresponding data set.

Please note that the functions used to compute the B-spline is taken from Milad Kharratzadeh’s post on B-splines. I have also tested the function various times, and I think it works properly, so that the problem is not with the definition of the B-spline.

``````functions {
//from: https://mc-stan.org/users/documentation/case-studies/splines_in_stan.html
vector build_b_spline(real[] t, real[] ext_knots, int ind, int order);
vector build_b_spline(real[] t, real[] ext_knots, int ind, int order) {
// INPUTS:
//    t:          the points at which the b_spline is calculated
//    ext_knots:  the set of extended knots
//    ind:        the index of the b_spline
//    order:      the order of the b-spline
vector[size(t)] b_spline;
vector[size(t)] w1 = rep_vector(0, size(t));
vector[size(t)] w2 = rep_vector(0, size(t));
if (order==1)
for (i in 1:size(t)) // B-splines of order 1 are piece-wise constant
b_spline[i] = (ext_knots[ind] <= t[i]) && (t[i] < ext_knots[ind+1]);
else {
if (ext_knots[ind] != ext_knots[ind+order-1])
w1 = (to_vector(t) - rep_vector(ext_knots[ind], size(t))) /
(ext_knots[ind+order-1] - ext_knots[ind]);
if (ext_knots[ind+1] != ext_knots[ind+order])
w2 = 1 - (to_vector(t) - rep_vector(ext_knots[ind+1], size(t))) /
(ext_knots[ind+order] - ext_knots[ind+1]);
// Calculating the B-spline recursively as linear interpolation of two lower-order splines
b_spline = w1 .* build_b_spline(t, ext_knots, ind, order-1) +
w2 .* build_b_spline(t, ext_knots, ind+1, order-1);
}
return b_spline;
}

//define exp(coeff*bspline)
//for testing purposes: just a single coeff/ b-spline
real exp_spline(real x, real xc, real[] theta, real[] x_r, int[] x_i)
{   real first;
first = x;
return exp(theta*build_b_spline(first, x_r, x_i, x_i));
}

//define a custom pdf - using the integral of an exponential b-spline
real espl_lpdf(real y, real theta, real coeff, data real[] ext_knots) {
real first;
real loc;
real deriv_spline;
real llik;
int x_i;
real theta_h;
first = y;
x_i = 1;
x_i = 2;
theta_h = coeff;
deriv_spline = build_b_spline(first, ext_knots, 1, 2);
print(ext_knots, y, theta_h, ext_knots, x_i);
loc = integrate_1d(exp_spline, ext_knots,y, theta_h, ext_knots, x_i, 1e-4);
llik = std_normal_lpdf(loc-theta)+coeff*deriv_spline;
return llik;
}
}

data {
int<lower=0> K;         // number of items
int<lower=1> N;         // number of test takers
matrix[N, K] Y;         // responses of test takers to the items
int num_knots;
real ext_knots[num_knots];
}

parameters {
vector[N] z;
vector[K] intercepts;
vector[K] coeffs;
real<lower=0> sigma_theta;
}

transformed parameters {
vector[N] theta = sigma_theta*z;
}

model {
sigma_theta ~ cauchy(0,1);
intercepts ~ normal(0,1);
coeffs ~ normal(0,1);
z ~ normal(0, 1);
for(i in 1:N)
{
for(j in 1:K)
{
Y[i,j] ~ espl(theta[i]-intercepts[j],coeffs[j],ext_knots); //
}
}

}
``````

Data file:

test_data.csv (5.4 KB)

In R, I then execute the following statements:

``````test_data <- read.csv("test_data.csv")
x <- as.vector(test_data)
minmax <- range(x)
irt_dat <- list(K = 3, Y=test_data, N=100, num_knots=4,ext_knots=c(minmax,minmax,minmax,minmax))
fit_spl <- stan(file = './stan_models/exp_spl.stan', data = irt_dat, iter = 4000, chains=1)
``````

This appears to be a bug in Stan 2.26 (and possibly earlier) which I can reproduce on my machine, but that does not appear with 2.27 where the model happily starts churning (and not really working unfortunately). Can you upgrade to the latest version and try again (latest rstan in unfortunately not on CRAN yet, but you can use Repository for distributing (some) stan-dev R packages | r-packages )?

Previous investigation if it turns out version is not a problem for future reference

So this is interesting. I modified the code to print all the points where the function gets evaluated and see what the function looks like when an error is thrown.

``````  real exp_spline(real x, real xc, real[] theta, real[] x_r, int[] x_i)
{   real first;
first = x;
print("   ",x," ",build_b_spline(first, x_r, x_i, x_i));
return exp(theta*build_b_spline(first, x_r, x_i, x_i));
}

``````

Here is one of the offending functions:

Code & data for the plot
``````to_plot <- tribble(~x, ~y,
-1.10844,0.596928
,1.846,0.213458
,-4.06289,0.980398
,1.99696,0.193865
,-4.21384,0.999991
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,0.985487,0.325148
,-3.20237,0.868708
,1.98931,0.194857
,-4.2062,0.998998
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,0.0629705,0.444885
,-2.27986,0.748971
,1.56092,0.25046
,-3.77781,0.943396
,1.95678,0.199079
,-4.17367,0.994776
,1.99606,0.193981
,-4.21295,0.999874
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,-1.10844,0.596928
,1.846,0.213458
,-4.06289,0.980398
,1.99696,0.193865
,-4.21384,0.999991
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,0.985487,0.325148
,-3.20237,0.868708
,1.98931,0.194857
,-4.2062,0.998998
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,0.0629705,0.444885
,-2.27986,0.748971
,1.56092,0.25046
,-3.77781,0.943396
,1.95678,0.199079
,-4.17367,0.994776
,1.99606,0.193981
,-4.21295,0.999874
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
,1.99703,0.193856
,-4.21391,1
)

to_plot %>%
ggplot(aes(x = x, y = exp(y * -0.367486))) + geom_line()
``````

Since Stan does under the hood integrate separately from lower bound to 0 and from 0 to upper bound, I also tried rewriting the function to shift the whole integration interval above zero, but this does not seem to have any effect.

1 Like

Thank you very much for your help. I did not yet change the version, but I am glad to hear that it resolves some bugs and I will try it. Thanks also for your efforts with the plotting and verification!

With respect to the plotted function: I deliberately started with a trivial toy-example. That is, the B-Spline I am using is indeed just a linear function. However, the exp() should make it nonlinear - unless the coefficient is very small (exp(coeff*bspline(x))), in that case the exp() might itself appear linear.
However, I really can’t understand why the quadrature has difficulties computing the integral of this (almost) linear function.