Question regarding integrate_1d example in docs

Modeling
1

Hi. In the docs, the example for where xc is useful in a 1D integral with finite limits, the model becomes:

real beta(real x, real xc, real[] theta, real[] x_r, int[] x_i) {
  real alpha = theta[1];
  real beta = theta[2];
  real v;
  if(x > 0.5) {
    v = x^(alpha - 1.0) * xc^(beta - 1.0);
  } else {
    v = x^(alpha - 1.0) * (1.0 - x)^(beta - 1.0);
  }
  return v;
}

But won’t there also be numerical instability near x=0 when \alpha is small?

Am I correct in assuming the following is valid?

real beta(real x, real xc, real[] theta, real[] x_r, int[] x_i) {
  real alpha = theta[1];
  real beta = theta[2];
  real v;
  if(x > 0.5) {
    v = x^(alpha - 1.0) * xc^(beta - 1.0);
  } else {
    v = xc^(alpha - 1.0) * (1.0 - x)^(beta - 1.0);
  }
  return v;
}

Is this a typo in the docs, or am I missing something blatantly obvious?

Thanks

2

xc is conceptually 1 - x. See

https://www.boost.org/doc/libs/1_71_0/libs/math/doc/html/math_toolkit/double_exponential/de_tanh_sinh_2_arg.html

2 Likes
3

The integral is over x and goes all the way to zero. The quadrature will evaluate points veeeery close to zero to get the integral.

If the sampler got to values of alpha that small, you might also have problems, but that isn’t guaranteed in the same way the integral evaluations will be (like, your posterior could be a long way away from zero, or more likely, the sampler never lands within like 1e-17 of alpha = 0.0 or whatever small number would trigger underflow).

1 Like
4

@mathDR or someone else, can you point me on that doc about xc? I can’t find it. Thanks!

5

Did you look at pg 191 of the manual? https://mc-stan.org/docs/2_22/stan-users-guide-2_22.pdf

2 Likes
6

Thanks - I just found it also (html: https://mc-stan.org/docs/2_22/stan-users-guide/integrator-convergence.html)

1 Like
7

Just bumping this. I think xc could also be 0 if x < 0.5 in the example. Is that correct?

8

This is the Boost docs which describes the behavior of the xc argument and what it’s good for: https://www.boost.org/doc/libs/1_73_0/libs/math/doc/html/math_toolkit/double_exponential/de_tanh_sinh_2_arg.html

We’re using their implementation directly so that info transfers.

1 Like
9

Yes, from that documentation, I infer that, if both endpoints are finite (assume [a, b]), then for arguments closer to a: xc = x-a, and arguments closer to b: xc = b-x.

So the Stan docs example could be rewritten as:

real beta(real x, real xc, real[] theta, real[] x_r, int[] x_i) {
  real alpha = theta[1];
  real beta = theta[2];
  real v;
  if(x > 0.5) {
    v = x^(alpha - 1.0) * xc^(beta - 1.0);
  } else {
    v = xc^(alpha - 1.0) * (1.0 - x)^(beta - 1.0);
  }
  return v;
}

Although I don’t (yet) know what the behavior is of xc when the distance to the endpoints is equivalent. (There must be a tie-breaking logic in there…)

1 Like
10

Yeah the model you wrote looks right. I think for a = 0 the complement logic isn’t necessary. x == xc for x < 0.5.

Regarding the switching point, I think that’s pretty arbitrary. The precision difference for this integral only matters close to the endpoints. You could probably change the logic to if(x > 0.01) and it wouldn’t matter.

1 Like
11

Well, nevermind lol. That’s definitely wrong cause the definition of xc is controlled by Boost. But I think Boost could move that threshold and it wouldn’t matter.