diff --git a/doc/NEWS.Rd b/doc/NEWS.Rd index ff7b983f78..4baa54dbd3 100644 --- a/doc/NEWS.Rd +++ b/doc/NEWS.Rd @@ -535,8 +535,8 @@ now, ditto for \code{arima0()}, thanks to \I{Norbert Kuder}'s report on the R-devel list. - \item \code{binomial()$linkinv(eta)} now also works for - \code{"logit"} link when \code{is.integer(eta)}. + \item \code{binomial()$linkinv(eta)} and \code{.. $mu.eta(eta)} + now also work for \code{"logit"} link when \code{is.integer(eta)}. } } } diff --git a/src/library/stats/src/family.c b/src/library/stats/src/family.c index 0e0399208d..542f0f3c9f 100644 --- a/src/library/stats/src/family.c +++ b/src/library/stats/src/family.c @@ -91,20 +91,18 @@ SEXP logit_linkinv(SEXP eta) SEXP logit_mu_eta(SEXP eta) { - int i, n = LENGTH(eta); - if (!n || !isReal(eta)) + int i, n = LENGTH(eta), nprot = 1; + if (!n || !isNumeric(eta)) error(_("Argument %s must be a nonempty numeric vector"), "eta"); + if (!isReal(eta)) {eta = PROTECT(coerceVector(eta, REALSXP)); nprot++;} SEXP ans = PROTECT(shallow_duplicate(eta)); double *rans = REAL(ans), *reta = REAL(eta); for (i = 0; i < n; i++) { - double etai = reta[i]; - double opexp = 1 + exp(etai); - - rans[i] = (etai > THRESH || etai < MTHRESH) ? DBL_EPSILON : - exp(etai)/(opexp * opexp); + double etai = reta[i], expE = exp(etai), opexp = 1 + expE; + rans[i] = (etai > THRESH || etai < MTHRESH) ? DBL_EPSILON : expE/(opexp * opexp); } - UNPROTECT(1); + UNPROTECT(nprot); return ans; } diff --git a/tests/reg-tests-1e.R b/tests/reg-tests-1e.R index a850b7f7c7..99177eec0b 100644 --- a/tests/reg-tests-1e.R +++ b/tests/reg-tests-1e.R @@ -1722,11 +1722,14 @@ stopifnot(exprs = { ## gave solve.default() error (as wrong model failed fitting) -## binomial()$linkinv() +## binomial()$ linkinv() and binomial()$ mu.eta() lnks <- c("logit", "probit", "cloglog", "cauchit", "log") binIlink <- function(eta) sapply(lnks, function(lnk) binomial(lnk)$linkinv(eta)) +binImuEt <- function(eta) sapply(lnks, function(lnk) binomial(lnk)$mu.eta (eta)) stopifnot(identical(binIlink( 0:3), binIlink(as.double(0:3)))) +stopifnot(identical(binImuEt( 0:3), + binImuEt(as.double(0:3)))) ## integer type was not allowed for logit (only) in R <= 4.4.2