# Using 1D numerical integration for doubly-intractable distributions

I am trying to fit the following four-parameter distribution in Stan, which is similar to the Subbotin distribution, but doesn’t have a closed-form expression. The kernel is:

f(x|\mu, T, \gamma, \phi) \propto \left(\frac{1}{1+e^{-\frac{x-\mu}{T}}}^{-\frac{1}{1+e^{-\frac{x-\mu}{T}}}}\frac{1}{1+e^{\frac{x-\mu}{T}}}^{-\frac{1}{1+e^{\frac{x-\mu}{T}}}}\right) e^{-\gamma \tanh \left(\frac{x-\mu}{2T}\right) (x-\mu)-\phi x}

The model is “doubly-intractable” in the sense that the normalizing constant, which is a function of the parameters, can’t be derived. I can get reasonable results writing my own Metropolis-Hastings algorithm that uses numerical integration to calculate the likelihood, and was wondering if it was possible to implement the integrate_1d function in Stan to do the same thing? I tried the following code:

functions{
real integrand(real x, real xc, real[] theta, real[] x_r, int[] x_i) {
real t = theta[1];
real mu = theta[2];
real gamma = theta[3];
real phi = theta[4];

return exp(-(1/(1+exp(-(x-mu)/t))*log(1/(1+exp(-(x-mu)/t)))+1/(1+exp((x-mu)/t))*log(1/(1+exp((x-mu)/t)))-
gamma*tanh((x-mu)/(2*t))*x-phi*(x)));
}
real model_lpdf(real x, real mu, real t, real gamma, real phi, real[] x_r, int[] x_i){
return -(1/(1+exp(-(x-mu)/t))*log(1/(1+exp(-(x-mu)/t)))+1/(1+exp((x-mu)/t))*log(1/(1+exp((x-mu)/t)))-gamma*tanh((x-mu)/(2*t))*x-phi*(x))-
log(integrate_1d(integrand, negative_infinity(),positive_infinity(),{mu, t, gamma, phi},x_r, x_i, 0.01));
}
}
data{
int<lower=0> N;
real x[N];
}
transformed data {
real x_r[0];
int x_i[0];
}
parameters {
real mu;
real<lower=0> t;
real<lower=0> gamma;
real phi;
}
model {
mu ~ normal(0,100);
t ~ cauchy(0,5);
gamma ~ cauchy(0,5);
phi ~ normal(0,100);
for(n in 1:N){
target += model_lpdf(x[n]  | mu, t, gamma, phi, x_r, x_i);
}
}


stan
Chain 1: Unrecoverable error evaluating the log probability at the initial value.
Chain 1: 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.
sinh_sinh quadrature cannot handle singularities in the domain.
If you are sure your function has no singularities, please submit a bug against boost.math
(in 'model160fb69b22a62_4d08a45af7a28f60f653fbfd036d9263' at line 17)
(in 'model160fb69b22a62_4d08a45af7a28f60f653fbfd036d9263' at line 41)


I was able to draw samples when I made the limits of integration finite (-10,10), but the chains would not converge and typically had a rejection rate of 95-100%. I suspect I’m using the integrate_1d function incorrectly or that Stan cannot accommodate such models. Any suggestions?

2 Likes

I would use expose_stan_functions (https://cran.r-project.org/web/packages/rstan/rstan.pdf) to experiment on this outside of the Stan model. Maybe there’s a bug in it, or maybe there’s a different way to code it that’s more stable or whatever.

return -(1/(1+exp(-(x-mu)/t))*log(1/(1+exp(-(x-mu)/t)))+1/(1+exp((x-mu)/t))*log(1/(1+exp((x-mu)/t)))- gamma*tanh((x-mu)/(2*t))*x-phi*(x));


The normalizing constant should be the integral of the function itself, not the log. The model lpdf would be the log density minus the log of that normalizing constant though.

1 Like

Thanks, and of course, oversight on expressing the integrand as the log.

If there is an actual singularity in the interior of the interval, then you would want to call integrate_1d twice, once with an upper bound at the singularity and again with the lower bound at the singularity.

3 Likes