Skip to content

Instantly share code, notes, and snippets.

@CGMossa
Last active March 15, 2024 11:05
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save CGMossa/955d1dbdcf203736111e1fa2a9c98b08 to your computer and use it in GitHub Desktop.
Save CGMossa/955d1dbdcf203736111e1fa2a9c98b08 to your computer and use it in GitHub Desktop.
There is a performance difference between invoking plogis(x) vs. binomial()$loglink(x), according to https://twitter.com/noah_greifer/status/1768505576143659132.
static const double MTHRESH = -30.;
SEXP logit_linkinv(SEXP eta)
{
SEXP ans = PROTECT(shallow_duplicate(eta));
int i, n = LENGTH(eta);
double *rans = REAL(ans), *reta = REAL(eta);
if (!n || !isReal(eta))
error(_("Argument %s must be a nonempty numeric vector"), "eta");
for (i = 0; i < n; i++) {
double etai = reta[i], tmp;
tmp = (etai < MTHRESH) ? DBL_EPSILON :
((etai > THRESH) ? INVEPS : exp(etai));
rans[i] = x_d_opx(tmp);
}
UNPROTECT(1);
return ans;
}
double plogis(double x, double location, double scale,
int lower_tail, int log_p)
{
#ifdef IEEE_754
if (ISNAN(x) || ISNAN(location) || ISNAN(scale))
return x + location + scale;
#endif
if (scale <= 0.0) ML_WARN_return_NAN;
x = (x - location) / scale;
if (ISNAN(x)) ML_WARN_return_NAN;
R_P_bounds_Inf_01(x);
if(log_p) {
// log(1 / (1 + exp( +- x ))) = -log(1 + exp( +- x))
return -log1pexp(lower_tail ? -x : x);
} else {
return 1 / (1 + exp(lower_tail ? -x : x));
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment