I’m trying to make a model with an exponentially modified Gaussian, with a scale parameter (not rate, as in the built-in version) that can vary between positive, zero, or negative.
In order to do this, I tried to implement the erfcx function, which seems to be missing from stan. Based on the example code at https://stackoverflow.com/questions/39777360/accurate-computation-of-scaled-complementary-error-function-erfcx , I’ve written the following functions:
functions {
real fmaf(real x, real y, real z) {
return x * y + z;
}
real erfcx(real x) {
real a;
real d;
real e;
real m;
real p;
real q;
real r;
real s;
real t;
a = abs(x); // NaN-preserving absolute value computation
/* Compute q = (a-2)/(a+2) accurately. [0,INF) -> [-1,1] */
m = a - 2.0;
p = a + 2.0;
r = 1.0 / p;
q = m * r;
t = fmaf (q + 1.0, -2.0, a);
e = fmaf (q, -a, t);
q = fmaf (r, e, q);
/* Approximate (1+2*a)*exp(a*a)*erfc(a) as p(q)+1 for q in [-1,1]
p = 0x1.f10000p-15f; // 5.92470169e-5
p = fmaf (p, q, 0x1.521cc6p-13f); // 1.61224554e-4
p = fmaf (p, q, -0x1.6b4ffep-12f); // -3.46481771e-4
p = fmaf (p, q, -0x1.6e2a7cp-10f); // -1.39681227e-3
p = fmaf (p, q, 0x1.3c1d7ep-10f); // 1.20588380e-3
p = fmaf (p, q, 0x1.1cc236p-07f); // 8.69014394e-3
p = fmaf (p, q, -0x1.069940p-07f); // -8.01387429e-3
p = fmaf (p, q, -0x1.bc1b6cp-05f); // -5.42122945e-2
p = fmaf (p, q, 0x1.4ff8acp-03f); // 1.64048523e-1
p = fmaf (p, q, -0x1.54081ap-03f); // -1.66031078e-1
p = fmaf (p, q, -0x1.7bf5cep-04f); // -9.27637145e-2
p = fmaf (p, q, 0x1.1ba03ap-02f); // 2.76978403e-1*/
p = 5.92470169e-5;
p = fmaf (p, q, 1.61224554e-4);
p = fmaf (p, q, -3.46481771e-4);
p = fmaf (p, q, -1.39681227e-3);
p = fmaf (p, q, 1.20588380e-3);
p = fmaf (p, q, 8.69014394e-3);
p = fmaf (p, q, -8.01387429e-3);
p = fmaf (p, q, -5.42122945e-2);
p = fmaf (p, q, 1.64048523e-1);
p = fmaf (p, q, -1.66031078e-1);
p = fmaf (p, q, -9.27637145e-2);
p = fmaf (p, q, 2.76978403e-1);
/* 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 = fmaf (p, r, r); // q = (p+1)/(1+2*a)
t = q + q;
e = (p - q) + fmaf (t, -a, 1.0); // residual: (p+1)-q*(1+2*a)
r = fmaf (e, r, q);
//if (a > 0x1.fffffep127f) r = 0.0f; // 3.40282347e+38 // handle INF argument
/* Handle negative arguments: erfcx(x) = 2*exp(x*x) - erfcx(|x|) */
if (x < 0.0) {
s = x * x;
d = fmaf (x, x, -s);
e = exp (s);
r = e - r;
r = fmaf (e, d + d, r);
r = r + e;
}
return r;
}
}
But when I try to compile this model, I get “Compilation ERROR, function(s)/method(s) not created!”. That’s a pretty opaque error. Can anybody help me debug it?