Repository navigation
Fix Exponentially modified Gaussian overflows at lambda >= 40 #2803
Description
Description
According to Wikipedia the exponentially modified Gaussian can be made more precise by a reparameterization and using the scaled erfc.
I have implemented the reparameterization and the issue below.
Example
library(Rcpp)
cppFunction('double my_erfc( double x){
return erfc( x ) ;
}')
cppFunction('double my_erfcx (double x)
// from https://stackoverflow.com/questions/39777360/accurate-computation-of-scaled-complementary-error-function-erfcx
{
double a, d, e, m, p, q, r, s, t;
a = fmax (x, 0.0 - x); // NaN preserving absolute value computation
/* Compute q = (a-4)/(a+4) accurately. [0,INF) -> [-1,1] */
m = a - 4.0;
p = a + 4.0;
r = 1.0 / p;
q = m * r;
t = fma (q + 1.0, -4.0, a);
e = fma (q, -a, t);
q = fma (r, e, q);
/* Approximate (1+2*a)*exp(a*a)*erfc(a) as p(q)+1 for q in [-1,1] */
p = 0x1.edcad78fc8044p-31; // 8.9820305531190140e-10
p = fma (p, q, 0x1.b1548f14735d1p-30); // 1.5764464777959401e-09
p = fma (p, q, -0x1.a1ad2e6c4a7a8p-27); // -1.2155985739342269e-08
p = fma (p, q, -0x1.1985b48f08574p-26); // -1.6386753783877791e-08
p = fma (p, q, 0x1.c6a8093ac4f83p-24); // 1.0585794011876720e-07
p = fma (p, q, 0x1.31c2b2b44b731p-24); // 7.1190423171700940e-08
p = fma (p, q, -0x1.b87373facb29fp-21); // -8.2040389712752056e-07
p = fma (p, q, 0x1.3fef1358803b7p-22); // 2.9796165315625938e-07
p = fma (p, q, 0x1.7eec072bb0be3p-18); // 5.7059822144459833e-06
p = fma (p, q, -0x1.78a680a741c4ap-17); // -1.1225056665965572e-05
p = fma (p, q, -0x1.9951f39295cf4p-16); // -2.4397380523258482e-05
p = fma (p, q, 0x1.3be1255ce180bp-13); // 1.5062307184282616e-04
p = fma (p, q, -0x1.a1df71176b791p-13); // -1.9925728768782324e-04
p = fma (p, q, -0x1.8d4aaa0099bc8p-11); // -7.5777369791018515e-04
p = fma (p, q, 0x1.49c673066c831p-8); // 5.0319701025945277e-03
p = fma (p, q, -0x1.0962386ea02b7p-6); // -1.6197733983519948e-02
p = fma (p, q, 0x1.3079edf465cc3p-5); // 3.7167515521269866e-02
p = fma (p, q, -0x1.0fb06dfedc4ccp-4); // -6.6330365820039094e-02
p = fma (p, q, 0x1.7fee004e266dfp-4); // 9.3732834999538536e-02
p = fma (p, q, -0x1.9ddb23c3e14d2p-4); // -1.0103906603588378e-01
p = fma (p, q, 0x1.16ecefcfa4865p-4); // 6.8097054254651804e-02
p = fma (p, q, 0x1.f7f5df66fc349p-7); // 1.5379652102610957e-02
p = fma (p, q, -0x1.1df1ad154a27fp-3); // -1.3962111684056208e-01
p = fma (p, q, 0x1.dd2c8b74febf6p-3); // 2.3299511862555250e-01
/* Divide (1+p) by (1+2*a) ==> exp(a*a)*erfc(a) */
d = a + 0.5;
r = 1.0 / d;
r = r * 0.5;
q = fma (p, r, r); // q = (p+1)/(1+2*a)
t = q + q;
e = (p - q) + fma (t, -a, 1.0); // residual: (p+1)-q*(1+2*a)
r = fma (e, r, q);
/* Handle argument of infinity */
if (a > 0x1.fffffffffffffp1023) r = 0.0;
/* Handle negative arguments: erfcx(x) = 2*exp(x*x) - erfcx(|x|) */
if (x < 0.0) {
s = x * x;
d = fma (x, x, -s);
e = exp (s);
r = e - r;
r = fma (e, d + d, r);
r = r + e;
if (e > 0x1.fffffffffffffp1023) r = e; // avoid creating NaN
}
return r;
}')
model_path <- cmdstanr::write_stan_file("
functions {
real my_exp_mod_normal_lpdf(real y, real mu, real sigma, real lambda) {
return exp_mod_normal_lpdf(y | mu, sigma, lambda);
}
}
")
expose_cmdstanr_functions(model_path, expose_to_global_env = T)
r_exp_mod_normal_lpdf <- function(y, mu, sigma, lambda) {
tau <- 1 / lambda
y_std <- (y - mu) / sigma
z <- (1 / sqrt(2)) * (sigma / tau - y_std)
logh <- -log(sigma) - 0.5 * log(pi)
half_log2pi <- 0.5 * log(pi) - log(2)
if (z < 0) {
x <- logh + log(sigma) - log(tau) + half_log2pi
return (x + 0.5 * (sigma/tau)^2 - ((y - mu) / tau) +
log(my_erfc(1/sqrt(2) * (sigma/tau - y_std))))
} else if (z >= 0 && z <= 6.71e7) {
logh <- -log(sigma) - 0.5 * log(pi)
x <- logh - 0.5 * y_std^2 + log(sigma) - log( tau)
return( x + half_log2pi + log(my_erfcx((1 / sqrt(2)) * (sigma / tau - y_std))))
}
return (logh - 0.5 * y_std^2 - log1p(y_std * tau / sigma))
}
> my_exp_mod_normal_lpdf(1, 0, 1, 40)
[1] -Inf
> r_exp_mod_normal_lpdf(1, 0, 1, 40)
[1] -1.394277
Activity
spinkney commentedon Aug 14, 2022
Having the
erfcxfunction in Stan will also help speed up the Wiener diffusion distribution and aid in the calculation of the cdf. See Even faster and even more accurate first-passage time densities anddistributions for the Wiener diffusion model by Matthias Gondan, Steven P. Blurton, Miriam Kesselmeier.
Also, if we implement a
log_Phifunction then usingerfcxcan be expanded for negative values. Below is the calculation for large negative values oflog(Phi(x))in stan-math which overflows and the calculation usingerfcx:This also can be used to make a
log_erffunction asspinkney commentedon Aug 14, 2022
Here's
log_Phiexpanded. It is equal tostd_normal_lcdfin the case of positive x (allows values of x up to 38) but expands negative x. Also see #2470.spinkney commentedon Oct 6, 2026
Closed by #3405