Using 1D numerical integration for doubly-intractable distributions

Modeling
1

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);
}
}

But received the following error,

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
2

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
3

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

4

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