p*() as the integral
of d*(), q*() inverting p*(), and the analytic gradients and Hessians
against numDeriv. They run in under a second.Log-space chain, subnormal band (utils.h, gkw.cpp, bkw.cpp,
kkw.cpp): the underflow bridge fired only at an exact 0, so for
x^alpha in [5e-324, 2.2e-308) ll*, gr* and hs* were silently wrong.
llbkw(c(161.8, 2, 1.5, 1), c(.01, .3, .6, .9)) was 1522.87 against a true
1523.21, and the gradient was off by up to 16%. The bridge now covers the band.
Ordinary data is bit-identical.
EKw had no bridge (ekw.cpp): dekw(), llekw(), grekw() and
hsekw() returned -Inf, +Inf or NaN where the nested GKw is finite, and
optim() stopped on them. They now match dgkw(gamma = 1, delta = 0).
Tails flushed to 0 or 1 in the p, q and r functions of GKw, BKw, KKw,
EKw, Kw and Mc. pgkw(1e-9, 40, 2, 0.05, 0.5, 0.1) returned 0 (true
0.0164), and rgkw() drew exact zeros that llgkw() then rejected. The
same held in the deep upper tail on the log scale (log.p = TRUE, log
probabilities below about -718). RNG streams are otherwise unchanged.
Upper tail of qgkw() and qmc(): both now reflect above y = 1/2, as
qbkw() already did. qgkw(1e-26, 2, 3, 1.5, 0.5, 1.2, lower.tail = FALSE)
returned exactly 1.
pgkw() and pbkw() with log.p = TRUE returned 0 near 1 for a tiny
negative log-probability.
hsgkw() rebuilt in log space: it returned NaN where llgkw() and
grgkw() are finite (e.g. beta = 200). It is now finite there, matches
numDeriv to 5e-9, equals hskkw() at gamma = 1, and is about 2x faster.
Memory leak on caught warnings: a warning raised from C++ and caught by
tryCatch() or options(warn = 2) skipped the C++ destructors (308 MB over
20 calls on 2e6 values). Warnings now unwind cleanly; messages are unchanged.
gkwgetstartvalues(family = "beta") started from Beta(gamma, delta + 2)
instead of Beta(gamma, delta + 1). Starting values are now the same on
every compiler.
Missing data give the documented value in all seven families: +Inf from
ll*(), NaN from gr*() and hs*().
safe_exp() no longer returns Inf for results between DBL_MAX / 10 and
DBL_MAX.
Analytic gradients and Hessians against numDeriv over 265 cases in all seven
families: maximum relative error 6.4e-10 and 5.4e-9 (before: 2e-2 and 1.0).
Two tests that pinned the old, wrong values were corrected.
gkw_log1mexp() had the sign of its second-order correction inverted
(utils.h): the Taylor branch, used for -1e-14 < u <= 0, returned
log(-u) - u/2 where the expansion gives log(-u) + u/2. Since
1 - exp(u) = -u (1 + u/2 + u^2/6 + ...), the correction is
log(1 + u/2) ~ +u/2, and the derivation in the comment above the line
carried the same slip. With u < 0 the two forms differ by |u|. Against a
400-digit reference:
u reference err before err after
-9.9e-15 -32.246241637770147 1.42e-14 0
-5.0e-15 -32.929338482476588 0 0
About 1.4 ulp at this magnitude, but in the wrong direction, and it left a
step where the function crosses into the expm1 branch -- breaking the
"relative error < 2*EPSILON" the header promises.
Downstream the effect is at the noise floor, and is reported as such: 56 of 238,140 grid values move, by at most 2.3e-13 relative, and adjudicating the changed densities against a 400-digit reference gives 4 closer and 2 farther, all between 1e-15 and 2e-14. The defensible claim is the direct one, on the function itself.
safe_pow() cast a double exponent to int without a range check
(utils.h): undefined behaviour once |y| > INT_MAX, which UBSan flags.
vec_safe_pow() had a guard for it, but that guard answered the parity
question wrongly in the process, reporting "even" for every |y| above
INT_MAX -- and an odd integer between INT_MAX and 2^53 is exactly
representable as a double. Both now take the parity from
fmod(|y|, 2) == 1, which is correct at every magnitude: above 2^53 every
double is even, and fmod says so.
safe_pow()'s documentation claimed an accuracy it does not have
(utils.h): it said the exp(y * log(x)) form "provides better numerical
stability than direct pow()". The reverse is true. That form carries a
relative error of roughly |y log x| * EPSILON, while std::pow on a
conforming libm is very nearly correctly rounded. Measured against a 60-digit
reference:
x y exp(y*log x) std::pow
10 100 1.11e-14 0
10 300 9.00e-14 0
2 1000 6.85e-14 0
The note now says what the form is actually for -- intercepting overflow and
underflow before they happen -- and points callers who do not need that to
std::pow.
The @return contract said NaN, the wrappers said stop() (all 28
d*, p*, q* and r* functions): the audit asked for one contract chosen
and propagated. stop() is the one kept -- it is what the package has shipped
for its whole CRAN life, a named error is a better diagnostic than a NaN for
the case that actually occurs, and switching to NaN would break every
tryCatch(..., error = ) guard written against it. All 28 @return blocks now
describe the behaviour the functions have. gr*() and hs*() genuinely do
return NaN, and their 14 blocks are unchanged.
Documentation only. No guard, message or value changes.
Known and deliberately left open: an infinite parameter is not intercepted
by the wrapper -- any(alpha <= 0) is FALSE for Inf and anyNA() does not
catch it -- so it reaches C++, where check_*_pars() has always rejected it,
and the routine leaves its fill value. Tightening the R guard is correct in
principle but is a behaviour change with a reverse-dependency cost:
gkwreg::predict() clamps only the lower bound of its linear predictor, so
exp(eta) can overflow to Inf on extrapolated newdata, and one such row
would abort the whole prediction vector. The value it gets today is 0, which is
the correct limit -- Kw(alpha, beta) concentrates at 1 as alpha grows -- so
nothing is currently wrong, only silent. An upper clamp in gkwreg is the
prerequisite. The current behaviour is pinned in
tests/testthat/test-return-contract.R so the eventual change is visible.
Only two of seven ll*() warned about data outside the open support
(llekw and llbeta): ll*() returns +Inf there -- an infinite objective
with no gradient direction, which an optimiser can only sit on -- and five of
the seven announced a corrupted sample not at all. All seven warn now.
Rcpp::warning() fired once per element inside the vectorised loops
(dgkw, pgkw, qgkw, rgkw, and the six nested r*): on 50,000 values
with half the parameters invalid that cost 51x, and under
options(warn = 2) each call longjmps out of the loop through C++ frames
holding live Armadillo objects. R's own convention is a single warning per
call. Each routine now sets a flag and warns once after the loop.
50,000 values, half invalid before after
elapsed 0.514 s 0.005 s
warnings raised 25,000 1
The message now also names what the routine actually returns, which base R
pairs up: dbeta(0.5, -1, 1) is NaN and warns "NaNs produced", while
rbeta(2, -1, 1) is NA and warns "NAs produced". p*, q* and r* fill
NA_REAL and say so. dgkw() leaves its fill value, 0, so it says
"invalid parameters" and claims no return value at all -- promising a NaN
that is not there is the defect this release fixed for q*().
Values are bit-identical over the 238,140-value regression grid. The path is
reachable from the exported API only through an Inf parameter, since every
wrapper stop()s on <= 0, NA and NaN.
Known inconsistency: man/dgkw.Rd said the function returns NaN for
invalid parameters; it returns 0, and the other six families document and
return 0 too. The @return contract entry above corrects the documentation.
d*() returned 0 at x = 0 and x = 1, where base R returns the limit
(all seven families): the GKw support is the open interval, but base R's
density functions carry the limiting value at the closed boundary --
dbeta(0, 0.5, 1) is Inf and dbeta(1, 2, 1) is 2 -- and any code that
plots a density across [0, 1] depends on it. Curves fell to zero exactly
where they should have diverged.
Substituting the first-order forms the log chain already uses, the log-density at each boundary collapses to a constant plus one power of the vanishing quantity, so a single exponent decides the answer:
at x = 0: alpha*gamma*lambda - 1 at x = 1: beta*(delta + 1) - 1
> 0 -> 0 = 0 -> the constant < 0 -> Inf
Each nested family reaches this through its own fixed parameters. The rule was
checked against stats::dbeta at all ten combinations the Beta
parameterisation can express, and against the nesting identities at both
boundaries for the other six.
dkw(0, 0.5, 1) 0 -> Inf dbeta_(1, 2, 0) 0 -> 2
Anything strictly outside [0, 1] is still 0, as in base R, and the interior
is bit-identical over the 238,140-value regression grid.
d*(), p*() and q*() dropped dim, dimnames and names (all seven
families): base R carries the first argument's attributes through to the
output, so code that indexes or plots the result by shape keeps working.
before after stats::dbeta
dim(d(matrix(x, 2, 2), 2, 3)) NULL 2 2 2 2
names(d(c(a = .2, b = .5), 2, 3)) NULL a b a b
The copy is conditional on the lengths agreeing, which is also what base R
does: once a recycled parameter makes the output longer than the first
argument, dim() is NULL in both. Values are bit-identical over the
238,140-value regression grid.
r*(0) raised an error instead of returning numeric(0) (all seven
families): stats::rbeta(0, 2, 3) is numeric(0), and a generator that
errors instead breaks any loop or replicate() that reaches an empty case.
All seven now return numeric(0). A negative or missing n is still an
error, as in base R.
pmc() reached R::pbeta through exp(lambda * log(x)), which is not a
round trip (bpmc.cpp): for lambda = 1 the Mc family is the Beta family,
and pmc(x, gamma, delta, 1) should be R::pbeta(x, gamma, delta+1) exactly.
It was not: exp(1 * log(x)) fails to return x for 2 of 9 ordinary values
under glibc, and for a different pair under the macOS ARM64 libm, which is
where CI caught it.
std::pow is used instead. C99 requires pow(x, 1.0) == x, so the identity
now holds bit-for-bit on every platform, and the form is also the more
accurate one for every other lambda. Against a 60-digit reference over 1,030
grid points:
lambda exp(l*log x) std::pow improved worse
1.5 1.19e-14 1.67e-16 403 3
0.1 6.57e-16 1.64e-16 56 0
2.5 1.57e-14 1.11e-16 527 0
Only pmc moves; 1,093 of 238,140 grid values, all in p mc.
The upper tail of pgkw() and pbkw() collapsed to exactly zero
(gkw.cpp, bkw.cpp): the CDF batch moved lower.tail and log.p onto
R::pbeta instead of applying them afterwards, which fixed the lower tail. A
second defect remained, in the argument rather than the tail flag.
pgkw() forms y = [1 - (1 - x^alpha)^beta]^lambda and evaluates
I_y(gamma, delta+1). As x approaches 1 the exponent lambda*log_w falls
below 1.1e-16, exp() returns exactly 1, and
R::pbeta(1, ., ., lower = FALSE) returns exactly 0. pbkw() reaches the
same place through -expm1 of an exponent running off to -Inf. Against a
300-digit incomplete beta, for pgkw(x, 2, 3, 1.5, 2, 0.8, lower.tail = FALSE):
1-x returned exact rel err
1e-4 5.731692e-34 5.731820e-34 2.2e-05
1e-5 5.840671e-43 5.734142e-43 1.9e-02
1e-6 0 5.734374e-52 1.00 <- collapse
1e-16 0 1.469535e-141 1.00
The true tail is still representable sixteen decades past the point where the
routine gave up, and lower.tail = FALSE is ordinary documented API, so
survival probabilities, p-values and quantile residuals were silently zero.
I_y(a,b) = 1 - I_{1-y}(b,a) is exact, and 1 - y is -expm1 of the same
exponent that produces y, at full relative accuracy. Reflecting above
y = 1/2 sends the small quantity into pbeta and the large one out of it;
below that crossover the direct form already holds the small quantity and is
left alone. This is the same correction pmc() received in this release.
After the fix the maximum relative error over those fifteen decades is
6.1e-14 for pgkw() and 5.3e-14 for pbkw(). Confinement over the
238,140-value grid: only pgkw and pbkw move, only with
lower.tail = FALSE, and no d* or q* value changes at all. The nesting
identity against pmc() -- corrected independently, in another translation
unit -- goes from 4.88e-01 to 2.84e-14.
NA, NaN and infinite input did not propagate (all seven families, 21
routines): NA_REAL is a NaN, so every comparison against it is false and
!R_finite() is true. Missing input therefore fell into the "outside the
support" branch and was silently replaced by the fill value:
before after base R
d*(NA) 0 NA NA
d*(NaN) 0 NaN NaN
p*(NA) 0 NA NA
p*(NaN) 0 NaN NaN
p*(+Inf) 0 1 1
p*(-Inf) 0 0 0
p*(Inf) = 0 also violated monotonicity outright: pkw(2, 2, 3) was 1 while
pkw(Inf, 2, 3) was 0. The fix needed no new branch -- dropping !R_finite()
from the lower boundary test lets +Inf fall through to the q >= 1 case
that was already there.
NA and NaN are distinguished, as base R distinguishes them: R_IsNA() is
asked before R_IsNaN(), so is.na() and is.nan() both answer correctly.
q*(p) saturated at 0 or 1 for probabilities outside [0, 1] (all seven
families): the R wrappers have always warned that such a value "will produce
NaN", while the C++ returned a bound. Defensive code testing is.nan() --
exactly what the warning tells the caller to expect -- saw nothing, and 0 and
1 are outside the open support, so the value flowed on into d*() and ll*()
as valid data.
qkw(-0.5, 2, 3) 0 -> NaN qkw(1.5, 2, 3) 1 -> NaN
qkw(-Inf, 2, 3) 0 -> NaN qkw(Inf, 2, 3) 1 -> NaN
The closed boundary is unchanged: q*(0) is still 0 and q*(1) still 1,
in both tails and on both scales.
Confinement was checked by bit-comparison rather than adjudication, since no in-range value should move: 238,140 values over the seven families, both tails, both scales, and x from 1e-320 to 1 - 1e-16 are bit-identical.
log, lower.tail and log.p silently accepted NA and read it as
TRUE (all seven families, 35 guards): the check was
!is.logical(flag) || length(flag) != 1, and NA passes both halves --
is.logical(NA) is TRUE and length(NA) is 1. The value reached C++ as
NA_LOGICAL, which is a non-zero integer, so it was read as TRUE:
dgkw(0.5, 2, 3, 1.5, 0.5, 2, log = NA) 0.8517722 <- the LOG density
dgkw(0.5, 2, 3, 1.5, 0.5, 2) 2.343797 <- the density
The documented error is now raised, with the message the help pages already promised.
An NA shape parameter behaved three different ways (all seven families,
92 guards): gkw, bkw, kkw, mc and kw wrote any(alpha <= 0), so if
received NA and R raised its own opaque "missing value where TRUE/FALSE
needed" instead of the documented message; ekw and beta_ wrote
any(alpha <= 0, na.rm = TRUE), which dropped the NA and returned 0 as
though the call had been valid. All seven now raise the package's own error:
dkw(0.5, NA, 3) before Error: missing value where TRUE/FALSE needed
dekw(0.5, NA, 3, 1) before 0
Both now raise the package's own error.
Valid calls are unaffected: dgkw(0.5, 2, 3, 1.5, 0.5, 2) is still 2.343797,
and delta = 0 is still accepted.
grkkw() and hskkw() failed silently and only partly (kkw.cpp):
alone among the seven families, they carried none of the guards their BKw
counterparts have. A refused input -- a short parameter vector, an invalid
parameter, data outside (0,1) -- returned a NaN result with nothing said,
so the caller could not tell a rejected call from a genuine boundary. Neither
function checked its result at all, so where the chain did reach a boundary
they returned whatever survived:
par = (1e-8, 1e-8, 0, 1e300), x = c(0.01, 0.3, 0.6, 0.9)
grkkw 1.378417e+307 -Inf -2.000000 31.43853 (no warning)
hskkw 2 of 16 entries NaN, 14 finite (no warning)
grbkw at the corresponding BKw point: warns, and every component is NaN
A gradient with a finite component next to an -Inf, or a Hessian with two
NaN entries among fourteen, is worse than no answer: it looks partly usable,
and the Hessian would have been inverted for a standard error. Both now warn
and return a uniformly NaN result, matching grbkw() and hsbkw().
The support test also gained has_nan(). A NaN compares false against both
bounds, so NaN data passed straight through; hskkw() returned a matrix
with fifteen NaN entries and one finite one.
This is a reporting change only. Over the same 43,792-value grid used above,
no value differs from the previous fix -- every case the guards catch was
already NaN or Inf.
The BKw and KKw likelihoods collapsed to +Inf on ordinary data
(bkw.cpp, kkw.cpp): both families walk the same log-space chain as their
GKw parent -- v = 1 - x^alpha, w = 1 - v^beta, z = 1 - w^lambda, with
lambda = 1 for BKw and gamma = 1 for KKw -- and both stopped at the point
where v underflows to exactly 1. log1mexp() then receives an argument of
exactly 0 and can only answer -Inf. The threshold is alpha*log(min(x)) < -745, which perfectly ordinary data crosses:
x = c(0.01, 0.3, 0.6, 0.9)
alpha = 161 llbkw(c(a,2,1.5,1), x) = 1515.5261017902 llkkw(c(a,2,1,1.5), x) = 1516.4186759096
alpha = 162 llbkw = Inf llkkw = Inf
alpha = 200 llbkw = Inf llkkw = Inf
Every larger alpha stayed at Inf, so the likelihood surface carried an
infinite plateau that an optimiser cannot leave. The true values, 1525.13 and
1526.03 at alpha = 162, are now returned.
The same boundary reached the density, the gradient and the Hessian.
dbkw() and dkkw() dropped such observations and returned the fill value
(-Inf in log, 0 otherwise). The derivatives built the ratios v^beta/w and
w^lambda/z as separate factors; each overflowed to +Inf on its own while
the quantity it multiplied had underflowed to 0, and 0 * Inf is NaN:
par = (200, 2, 1.5, 1), x = c(.01, .30, .60, .90)
grbkw NaN NaN NaN NaN
grgkw(lambda=1) 9.617994 -3.000000 1278.026571 -2.721489
par = (200, 2, 1, 1.5), x = c(.01, .30, .60, .90)
grkkw NaN NaN -2 NaN (partly NaN, and silently so)
grgkw(gamma=1) 9.617994 -3.000000 -2.000000 1279.626571
Both files now use the first-order limits that are exact to the last
representable bit -- as x -> 0, log_w = log(beta) + alpha*log(x); as
x -> 1, log_z = log(lambda) + beta*log_v -- keep every ratio inside a
single exp() of a sum of logs, and never let a coefficient of exactly zero
multiply a logarithm. This is the same repair gkw.cpp received, so the
nesting identities BKw(a,b,g,d) == GKw(a,b,g,d,1) and
KKw(a,b,d,l) == GKw(a,b,1,d,l) now hold where they previously did not.
Over a grid of 19,488 values spanning x from 1e-300 to 1 - 1e-16 and
alpha from 0.01 to 5000, 1,966 results changed from NaN/Inf to a finite
value and none went the other way. 1,089 NaN Hessian entries -- 71 whole
unusable matrices -- became finite. Adjudicated against a 400-digit reference
and against the GKw parent: the maximum relative error fell from infinite to
1.8e-07 for the log-likelihoods and 1.4e-04 for the gradients, both attained
only at alpha >= 1000 where double precision has no digits left to lose.
No value that was already exactly correct changed; values correct to within
one ulp rose from 100 to 152 (log-likelihood), 348 to 479 (gradient) and 1,337
to 1,716 (Hessian). pbkw(), qbkw(), pkkw() and qkkw() are untouched
and bit-identical.
pmc()'s upper tail was quantised by the argument it handed to
R::pbeta (bpmc.cpp): F(x) = I_{x^lambda}(gamma, delta+1), and pmc()
formed x^lambda in linear arithmetic. Once x^lambda passes 1/2 a double
holds it no more finely than 1.1e-16, and the upper tail is a function of
1 - x^lambda alone, so it was quantised to whatever that left -- and to
exactly 0 once x^lambda reached 1:
pmc(x, 4, 2, 0.25, lower.tail = FALSE) before exact
x = 1 - 1e-15 2.1895288505e-46 3.1175127579e-46
x = 1 - 1.1e-16 0 4.2764235361e-49
a relative error of 700%, then of 100%. I_y(a,b) = 1 - I_{1-y}(b,a) is
exact, and 1 - x^lambda comes from -expm1 of the same exponent at full
relative accuracy, so reflecting sends the small quantity into pbeta. The
reflection is applied only where the direct form is the one holding the large
quantity: the lower tail never changes, and neither does any upper tail with
x^lambda <= 1/2.
Over 8942 grid cells whose exact tail is representable at all, 8494 are
bit-identical, 353 improved and 76 moved the other way. The maximum relative
error fell from 7.00 to 1.36e-13, and that residual sits on a cell this change
did not touch. The 76 that moved the other way went from at most 2.67e-14 to
at most 4.10e-14 relative, which is R::pbeta's own accuracy at tail
probabilities of 1e-71: both routes hand it an argument good to one ulp there,
and they round differently. dmc(), qmc(), rmc(), llmc(), grmc() and
hsmc() are bit-identical, as is every lower tail and every lambda = 1
result, which stays identical to stats::pbeta to the bit.
grmc() and hsmc() swapped R::digamma and R::trigamma for two-term
asymptotic expansions above three separate thresholds (bpmc.cpp):
gamma > 100, delta > 100 and gamma + delta > 100. log(z) - 1/(2z)
truncates psi's expansion before the 1/(12z^2) term and is wrong by 8.33e-06
at z = 100 and by 1.30e-03 at z = 8, which the gamma + delta threshold
can reach with gamma that small; 1/z + 1/(2z^2) drops psi'-s 1/(6z^3)
term and is wrong by 1.67e-07 at z = 100.
Each threshold put a step of n times that error into a different component,
at a different place. On the seven observations c(.1,.25,.4,.5,.6,.75,.9):
grmc(c(gamma, 3, 1), x)[1] gamma = 99.999 5.9262333798
gamma = 100.001 5.9262971512
a jump of 6.38e-05, of which 5.83e-05 is discontinuity rather than slope; the
step scales with n and reaches 0.018 at n = 2160. The largest step between
neighbouring gamma falls from 6.11e-05 to 2.72e-06 in the gradient and from
1.22e-06 to 5.38e-08 in H[gamma, gamma]. R::digamma and R::trigamma are
accurate to 1e-16 over the whole range, so the substitution bought nothing.
Over 330 gradient cells the maximum relative error fell from 3.85e-04 to
1.07e-12 and no cell got worse; over 660 Hessian cells 105 improved and 20
moved the other way. Ten of those twenty are H[gamma, delta] moving by one
ulp. The other ten are H[gamma, gamma] and H[delta, delta] at
gamma or delta = 1e12, where the entry is a difference of two psi' values
that agree to eleven digits: the result, around 1.8e-23, can carry no better
than 1e-04 relative in double precision whatever psi' returns, and the
measured relative error moves from 1.9e-05 to 7.2e-04, both inside that floor.
The smooth asymptotic form landed inside it by luck at that one point while
being 5.1e-04 wrong and discontinuous at the far more plausible gamma = 100.
dmc(), pmc(), qmc(), rmc() and llmc() are bit-identical.
llmc() swapped R::lbeta for a difference of lgamma above
gamma = 100 or delta = 100 (bpmc.cpp): that difference is the
cancellation R::lbeta exists to avoid. At gamma = 1e12, delta = 2 the
two outer lgamma values are 2.66e13, where one ulp is 3.9e-03, so their
difference cannot resolve an answer of -82.2 any better than that:
R::lbeta(1e12, 3) -82.1999161672287 (exact)
lgamma(1e12) + lgamma(3) - lgamma(1e12+3) -82.203125 (off by 3.2e-03)
dmc() and llbeta() always called R::lbeta, so llmc() also disagreed
with -sum(dmc(..., log = TRUE)), the objective it is supposed to be, and
with llbeta() at lambda = 1, where the two are the same model:
llmc(c(1e12, 2, 1), 1 - 1e-12) before -25.9371543959 exact -25.9378518132
llmc(c(2, 1e12, 1), 1e-12) before -26.6306976341 exact -26.6310211159
R::lbeta is now called at every gamma and delta. Over 110
likelihood cells the maximum error fell from 16283 ulps of the working
magnitude to 3.00, which is where the branch-free regimes already sat; 91
cells are bit-identical, 16 improved and 2 moved by at most 2 ulps.
dmc(), pmc(), qmc(), rmc(), grmc() and hsmc() are bit-identical.
dmc() lost the density as x approached 1 (bpmc.cpp): it formed
x^lambda in linear arithmetic and then took log(1 - x^lambda). Doubles are
spaced 2.2e-16 apart just below 1, so 1 - x^lambda carries an absolute error
of one ulp of 1 however small it truly is, and x^lambda rounds to exactly 1
once 1 - x drops under about 1e-16, at which point a guard returned a density
of zero.
Mc(gamma, delta, lambda) is GKw(1, 1, gamma, delta, lambda), so dgkw()
with alpha = beta = 1 is the same density computed a different way, and it
was already right:
gamma = 1.5, delta = 2, lambda = 0.8
exact (400 digits) dgkw(x,1,1,..) dmc
x = 1 - 1e-13 -58.65464965 -58.65464965 -58.65409479
x = 1 - 1e-15 -67.86721101 -67.86721101 -67.92355276
x = 1 - 1e-16 -72.26166017 -72.26166017 -71.81537306
log(1 - x^lambda) now goes through gkw_log1mexp(lambda * log(x)), the
helper dgkw() already uses, and the guard is gone.
llmc() and grmc() carried the mirror image of the same defect at the other
end of the support: log(-expm1(u)) has to represent a number just below 1
and so reported log(1 - x^lambda) as a multiple of 1.11e-16 -- usually as
exactly 0 -- for every x^lambda under one ulp. At delta = 1e12 the missing
term is worth 2e-05 nats an observation. All four functions now share
gkw_log1mexp().
Adjudicated against a 120-digit decimal reference over 9130 density cells
(22 parameter settings x 415 quantiles from 5e-324 to 1 - 1.1e-16): 8424 cells
are bit-identical, 512 improved, 143 moved by at most 4 ulps of the working
magnitude, and the maximum error fell from infinite -- six cells returned
-Inf for a finite density -- to 10 ulps. The nesting identity
dmc(x, g, d, l) == dgkw(x, 1, 1, g, d, l) closed from 9.8e13 ulps to 10.
pmc(), qmc(), rmc(), hsmc(), dgkw(), llgkw() and llbeta() are
bit-identical across the whole grid.
dgkw(), llgkw() and grgkw() broke down along the same log-space
chain (gkw.cpp): all three walk v = 1 - x^alpha, w = 1 - v^beta,
z = 1 - w^lambda, and each lost the chain in its own way.
llgkw() computed log(x^alpha) as vec_safe_log(vec_safe_pow(x, alpha)), a
round trip that both lost digits and made it disagree with dgkw(), which
already used alpha * log(x).
Two of the three transformations underflow to a boundary that log1mexp()
cannot recover from: its argument arrives as exactly 0 and log(1 - exp(0))
is -Inf. With a zero coefficient in front -- delta = 0, or
gamma*lambda = 1 -- 0 * -Inf is NaN:
llgkw(c(1, 300, 1, 0, 1), c(.8, .85, .9, .95)) NaN, for an exact 2609.84
Both regimes have a first-order limit that is exact to the last representable
bit: as x -> 0, log_w = log(beta) + alpha*log(x); as x -> 1,
log_z = log(lambda) + beta*log_v. These are now used where the direct form
underflows, and a coefficient of exactly zero never multiplies a logarithm.
grgkw() built 1/v, 1/w and 1/z as separate reciprocals, each of which
overflowed on its own long before the product it belonged to was large:
par = (1, 70, 1.5, 2, 1), x = c(.10, .25, .40, .72, .99)
llgkw 1386.78983044 (finite, correct)
grgkw NaN NaN NaN NaN NaN
numDeriv::grad -662.27 20.27 -6.76 472.41 -15.00
Correcting LOG_DBL_MAX in 1.1.6 moved that boundary out by 2.3x but did not
remove it; beta = 70 recovered while beta = 200 still failed. Every ratio
is now a single exp() of a difference of logs, so only the difference has to
be representable, not the reciprocal.
Over 96 parameter/data blocks: llgkw() was non-finite in 30 and is now
non-finite in none; grgkw() returned NaN in 30 and now in none. Adjudicated
against a 900-digit reference and against grbkw(), an independent
implementation of the same gradient at lambda = 1: the maximum error in
llgkw() fell from 8.7e-05 to 7.5e-16, and grgkw() agrees with grbkw() to
5.4e-08. dgkw() dropped 812 spurious -Inf log-densities across eight
parameter settings while leaving every already-finite value bit-identical.
Note on references: numDeriv is not usable as an arbiter in part of this
region. For beta = 1000 the intermediate log_w becomes subnormal, with 15
significant bits left, and lambda * log_w is then bit-identical for
lambda = 1 +/- 1e-6 -- the finite difference sees no dependence at all and
reports a lambda component short by exactly delta/lambda. The analytic
value is correct there; grbkw() and the 900-digit reference confirm it.
grkw() and hskw() were not the gradient and Hessian of llkw()
(kw.cpp): kw.cpp was the last family file whose derivatives were still
evaluated in linear arithmetic. It formed v = 1 - x^alpha and then applied
arma::clamp(v, eps, 1 - eps) with eps = 2.22e-14, freezing log(v) at
-31.4384832 for every observation near 1 regardless of the data:
par = (0.5, 2), x = 1 - 1e-14 d/dalpha d/dbeta
grkw -2.44999999999999 30.938483203129
numDeriv::grad(llkw) -3.99999999995549 32.430138079912
grekw(lambda = 1), same dist. -3.99999999999997 32.430138079907
The relative error reached 38.75% in the gradient, 47% in H[alpha,alpha] and
78% in H[alpha,beta], and ordinary data was enough to show it:
c(1-1e-9, 1-1e-11, 0.5) already diverged by 3.4e-05. Standard errors and
confidence intervals were wrong whenever the sample held observations close
to 1.
Both routines are now grekw() / hsekw() with lambda fixed at 1, evaluated
from logarithms. Adjudicated over 80 parameter/data blocks against two
independent references: the maximum relative disagreement with the EKw path
fell from 3.93 to exactly 0 in the gradient and from 4.00 to 1.4e-14 in the
Hessian, and against numDeriv from 3.93 to 1.6e-09, which is numDeriv's own
noise. No block moved in the wrong direction and no sign changed.
llgkw() returned -Inf for data outside the open support (gkw.cpp):
ll*() is the negative log-likelihood, so an invalid point must be +Inf --
the value optim() moves away from. llgkw() returned -Inf, making data
outside (0, 1) the global minimum of the function being minimised, and it
was the only one of the seven families with that sign; llbkw(), llkkw(),
llekw(), llkw(), llmc() and llbeta() all returned +Inf. The parameter
path in llgkw() was already +Inf and is unchanged.
The practical damage was in comparing likelihoods rather than in optimisation:
optim() refuses to start at either infinity, so no optimiser silently
converged on bad data. But on a sample holding a single 0 -- ordinary in
untransformed proportions -- the GKw family won every comparison:
gkw nll = -Inf <- wins argmin = gkw
bkw nll = Inf AIC = -Inf
... nll = Inf
Values for valid data are bit-identical, and the nesting identities against
llkw(), llbkw() and llekw() still agree to 1e-12.
Zero-length arguments crashed the R process (all seven families): the
vectorised d*(), p*(), q*() and r*() routines size their output as the
maximum length of their inputs and then recycle with i % vec.n_elem. When one
argument had length zero while another did not, the output length stayed at one
or more and the recycling evaluated i % 0. Integer division by zero is
undefined behaviour; on x86-64 it raises SIGFPE, terminating the R process
with no error, no message and nothing for tryCatch() to catch. The R-level
validation did not intercept it either, because any(numeric(0) <= 0) is
FALSE. All 28 exported routines now short-circuit before the loop:
d*(), p*() and q*() return numeric(0), matching
stats::dbeta(numeric(0), 1, 1), and r*() return n missing values with a
warning, matching stats::rbeta(3, numeric(0), 1). A filtered vector that
happened to be empty, such as dkw(x[x > 1], 2, 3), was enough to trigger the
crash. Numerical output is unchanged for every non-empty input.
rbkw() and rkkw() generated values outside the open support
(gkw.cpp, bkw.cpp, kkw.cpp, ekw.cpp, kw.cpp): rbkw() drew
V ~ Beta(gamma, delta+1) and then formed 1.0 - V. For V below 1.1e-16
that rounds to exactly 1 and the generator returned 0. R::rbeta itself never
returned a zero -- every one was fabricated by the subtraction:
rbkw(1e5, 2, 3, 0.02, 0) 48,602 exact zeros 48.6% of the sample
rbkw(1e5, 2, 3, 0.05, 0) 16,440 exact zeros
rkkw(1e5, 0.2, 3, 0, 0.3) 6 exact zeros
A zero is outside the (0,1) the likelihood accepts, so the package's own
simulate-then-fit workflow broke on six zeros in a hundred thousand:
llkkw() at the true parameters returned Inf and optim() stopped with
"L-BFGS-B needs finite values of 'fn'". Kolmogorov-Smirnov against the
package's own CDF rejected the gamma = 0.02 sample outright, D = 0.486.
The five generators that formed 1 - u -- rgkw(), rbkw(), rkkw(),
rekw(), rkw() -- now invert in log space, the same chain the quantile
functions use. rmc() and rbeta_() never had the defect and are untouched.
The draws themselves are unchanged, so set.seed() reproduces exactly the
stream it did before: .Random.seed after 100,000 variates is bit-identical
for all seven generators. Only the inversion that follows differs. Replaying
the same Beta draws, 97,559 of 100,000 rbkw() values changed and every one
moved closer to the closed-form inversion, none away; the largest relative
error falls from 1.0 to 7.2e-15, and for rkkw() from 1.0 to 8.5e-14.
The Kolmogorov-Smirnov statistic for gamma = 0.02 goes from D = 0.486
(p < 1e-16, 24,277 variates outside the support) to D = 0.0042,
p = 0.336, none outside. The fit now converges and recovers the parameters.
All seven quantile functions returned values outside the open support
(gkw.cpp, bkw.cpp, kkw.cpp, ekw.cpp, bpmc.cpp, kw.cpp,
beta_.cpp): each q*() undid log.p with exp(), folded the upper tail
with 1 - p, and then inverted using 1 - u in linear space at every step.
The result was not merely imprecise -- it left (0,1) altogether:
qekw(0.02, 20, 0.1, 0.1) returned 0 true value 0.1587
qekw(1e-08, 5, 2, 0.5) returned 0 true value 5.49e-04
qkw(1e-16, 0.2, 2) returned 0 true value 3.125e-82
qbeta_(-1000, 2, 3, log.p=T) returned 0 true value 2.25e-218
A quantile of exactly 0 or 1 then feeds d*() and ll*() a value outside the
support they accept, so the damage propagates into simulation by inversion and
into any likelihood built on it.
The inversion now carries log(u) and log(1-u) from whatever scale and tail
the caller used, so neither is recovered by subtraction, and the four families
that route through the incomplete beta hand lower_tail/log_p to
R::qbeta. qbkw() needs log(1-z): it takes log1p(-z) while z <= 1/2
and otherwise gets 1-z directly from R::qbeta through the symmetry
I_z(a,b) = 1 - I_{1-z}(b,a).
Judged against closed-form inversions written independently in R: of 22,236
values, 8,311 changed and every one improved. The largest relative error falls
from 2.06 to 5.7e-14. The boundary conventions of 1.1.5 are preserved
exactly, including the saturating result for out-of-range p.
One limit remains: below about log(p) = -745, exp(log u) underflows, 1-u
rounds to exactly 1 and the inversion has nothing left to invert, so the
quantile is still 0. Recovering it needs each step to carry both log(q) and
log(1-q). 1.1.5 already returned 0 from log(p) = -40 downward.
All seven cumulative distribution functions collapsed to 0 or 1
(gkw.cpp, bkw.cpp, kkw.cpp, ekw.cpp, bpmc.cpp, kw.cpp,
beta_.cpp): each p*() formed 1 - x^alpha and 1 - (1 - x^alpha)^beta
in linear space. Once x^alpha fell below 1.1e-16 the first rounded to
exactly 1 and the second to exactly 0, and the CDF returned 0 or 1. The error
was absolute, not a lost digit:
pekw(5.62e-09, 2, 5, 0.02) returned 0 true value 0.483
pekw(0.14, 20, 20, 0.1) returned 0 true value 0.0264
pkw(1e-09, 2, 5, log.p=TRUE) returned -Inf true value -39.84
The second sits at x = 0.14, nowhere near a tail: a small lambda
compresses the result toward 1 and pulls the collapse into the body of the
distribution. A systematic sweep found 8,917 affected points for pgkw()
alone, with a maximum absolute error of 1.0 -- the largest a probability can
be wrong by.
lower.tail and log.p were also applied afterwards, as 1 - p and
log(p), instead of being passed to R::pbeta, which implements both without
ever forming those quantities. That cost the opposite tail:
pmc(1 - 1e-06, 2, 3, 2.5, lower.tail = FALSE) returned exactly 0 against a
true 1.95e-22, and pbeta_(1e-200, 2, 3, log.p = TRUE) returned -Inf
against a true -918.73.
Every chain now runs in log space through gkw_log1mexp(), the survival
function is computed directly rather than as 1 - F, and lower_tail/log_p
go straight to R::pbeta where the family routes through the incomplete beta.
Judged against log-space references written independently in R: of 61,100
values over a grid spanning 1e-300 to 1 - 1e-16, 23,720 changed, 23,719
improved and one moved by a single ulp. The largest relative error falls from
5.9e+305 to 1.6e-01, and only 14 values of the 61,100 still exceed 1e-9.
Those 14 sit at x >= 0.999, where x^lambda is within an ulp of 1 and the
argument handed to R::pbeta cannot carry more precision. All seven CDFs are
monotone, stay within [0, 1], satisfy F + S = 1 to 1.1e-16, and
reproduce their nesting identities exactly.
The McDonald log-likelihood was unbounded below, and silently rewrote the
data (src/bpmc.cpp): llmc(), grmc() and hsmc() clamped every
observation to [1e-10, 1-1e-10] before use. For a legitimate observation at
1e-20 that moved the likelihood by 23 nats, and it broke the identity
llmc(gamma, delta, 1) == llbeta(gamma, delta), where the two are the same
model: the disagreement reached 1140 nats.
Separately, for delta > 1000 the term delta * log(1 - x^lambda) was
floored at -700 per observation, so it stopped growing with delta while
the constant term n(log lambda - log B(gamma, delta+1)) kept growing. The
negative log-likelihood became unbounded below:
llmc(c(1e300, 1e300, 1e-6), x) returned -2.77e+302, a global minimum at
absurd parameters, with a visible step at delta = 1000. It now returns
+2.58e+303.
grmc() and hsmc() additionally floored v = 1 - x^lambda at 1e-10 and
capped their lambda terms at ±1e6, so the gradient plateaued where the
objective kept moving. llmc() also computed log(1 - x^lambda) as
log1p(-x^lambda), which cannot recover digits x^lambda has already lost,
while grmc() and hsmc() already used -expm1(lambda * log(x)); the
objective and its gradient therefore disagreed as x approached 1. All three
now use -expm1 of the same exponent.
These clamps were removed together rather than one at a time: llmc() shared
them with grmc() and hsmc(), so removing only the documented subset would
have introduced an objective/gradient mismatch that did not previously exist.
Because the three functions shared the same clamps, checking the analytic
gradient against numDeriv::grad(llmc) passed even with the defect present.
Validation therefore used a closed-form reference written independently in R.
Against it, the largest relative error falls from Inf to 1.0e-14 for
llmc(), from 1.12 to 3.1e-05 for grmc(), and hsmc() now agrees with
the jacobian of the analytic gradient to 6.0e-08. The nesting identity with
stats::dbeta() goes from 1140 nats of error to 9.1e-13. Maximum-likelihood
fits on well-behaved data are bit-identical, under both BFGS with the analytic
gradient and Nelder-Mead.
The residual 3.1e-05 in grmc() appears only for gamma + delta > 100 and
is unchanged by this commit: it comes from the asymptotic digamma expansions,
a separate defect.
hsgkw() silently returned the Hessian of a smaller sample
(src/gkw.cpp): the observation loop skipped past any point whose
log(1-x^alpha), log(1-v^beta) or log(1-w^lambda) came out non-finite,
leaving the remaining terms to be returned as a finite, symmetric matrix with
no NaN and no warning. For a quantity whose purpose is to produce standard
errors, that is the worst available failure mode. With beta = 500 and four
observations every point was dropped and only the parameter-only terms
survived, so H(alpha, alpha) came back as n / alpha^2 = 4 against a true
1996.3 -- wrong by a factor of 499, and indistinguishable from a valid
result. With five observations one survived and the function returned the
Hessian of a single point as if it described all five.
The loop now stops on the first such observation and returns a NaN matrix
with a warning, matching what the function's own intermediate-value check
already did. Matrices are bit-identical wherever they were finite before.
Computing those terms correctly requires the log-space rework and is not
attempted here. grgkw() is unaffected: it is fully vectorised and propagates
NaN rather than dropping observations.
dgkw() returned a density of zero as x approached 1 (src/gkw.cpp):
the density formed x^alpha in linear space and bailed out whenever
x^alpha >= 1 - sqrt(.Machine$double.eps). The guard was there because
log(x^alpha) loses its significant digits in that band -- doubles are spaced
2.2e-16 apart near 1, so the relative error reaches 4e-6 by
1 - x = 1e-12 -- but returning zero is a far worse answer than an imprecise
one. It also broke the nesting identity: dgkw(1 - 1e-9, 1, 0.1, 1, 0, 1)
returned 0 while dkw(1 - 1e-9, 1, 0.1), the same density, returned
1.26e+07. For GKw(0.1, 0.1, 10, 0.1, 0.1) the discarded band held 13% of
the probability mass, and for beta < 1, where the density diverges at 1, the
rising tail was replaced by a cliff to zero.
log(x^alpha) is now taken as alpha * log(x), which is exact and removes the
need for the guard; gkw_log1mexp() already covers the resulting regime. Over
the regression grid, 6,780 of 47,104 dgkw() values changed, every one closer
to an independent log-space reference and none further away; no other family
moved. Recovered mass shows up in the integral: GKw(0.1, 0.1, 10, 0.1, 0.1)
goes from 0.8671 to 0.9998.
Two continue guards further down dgkw() still discard a point when
log(1 - w^lambda) underflows, even where delta = 0 makes that term vanish
from the density. That is unchanged here and belongs with the log-space rework.
Two logarithmic bound constants held the wrong quantity (src/utils.h):
LOG_DBL_MAX was documented as log(DBL_MAX_SAFE) but held
log10(DBL_MAX) = 308.2547, while the correct natural logarithm is 707.4801.
Since it is used as the overflow threshold of safe_exp() and safe_pow(),
every result above exp(308.25) was returned as +Inf, discarding roughly 174
orders of magnitude of representable double range. Reachable from the public
API: dkw(1e-300, 0.5, 2) returned Inf instead of 1e+150.
Separately, safe_log() scaled its underflow branch by LOG_DBL_MIN, which is
log(DBL_MIN), while dividing by DBL_MIN_SAFE, which is 10 * DBL_MIN. Every
result below 2.225e-307 was therefore off by exactly log(10) = 2.302585 --
a finite, plausible, wrong number rather than a visible failure. It propagated
into dkw(x, log = TRUE), llkw(), llgkw() and pmc(log.p = TRUE), and made
llgkw() disagree with dgkw(), which takes log(x) directly.
The constants are now named for what they are -- LOG_DBL_MIN,
LOG_DBL_MIN_SAFE and LOG_DBL_MAX -- and safe_log() scales by the logarithm
of the divisor it actually uses. Over a regression grid of 401,373 values, 328
density values, 16 log-likelihoods and 3 tail probabilities changed; every one
moved closer to an independent log-space reference, and none moved away. The
largest relative error against that reference fell from Inf to 1.6e-16 for
densities and from 4.5% to 1.2e-10 for llgkw(). As a side effect the
all-NaN region of grgkw() recedes: with x_max = 0.99 it began at
beta = 80 and now extends past beta = 130, with the newly finite values
agreeing with numDeriv::grad() to 1.3e-9 or better.
safe_exp() still saturates above log(DBL_MAX_SAFE), i.e. one order of
magnitude below the true double maximum. That headroom is the documented intent
of the DBL_MAX_SAFE constant and is left in place.
gkwgetstartvalues() never ran the multi-start it documents; n_starts was
inert (gkwinit.cpp): the selection loop decided whether to optimize a
starting point by comparing the raw objective at that point against
best_obj, which after the first iteration already held an optimized value.
A start that was merely poor -- exactly the case multi-start exists to rescue --
was therefore discarded before Nelder-Mead ever saw it, and in practice only
the first of the candidates was optimized at all. Every value of n_starts
returned the same answer, bit for bit:
set.seed(202); x <- rgkw(500, 5, 1.2, 3, 0.5, 2)
gkwgetstartvalues(x, "gkw", n_starts = k) moment objective
k = 1 6.262579e-04
k = 10 6.262579e-04 (identical)
k = 200 6.262579e-04 (identical)
after
k = 1 1.008584e-07
k = 10 1.369536e-08
k = 200 1.037278e-10
Every candidate is now optimized and only the optimized objectives are
compared. A second defect surfaced once the multi-start was live: Nelder-Mead
is unconstrained and can leave the family's parameter box, so a winner chosen
on its pre-clamp objective could be handed back worse than a rival that stayed
inside, making the returned error non-monotone in n_starts. Candidates are
now clipped to the box before being scored, so the objective that is compared
is the objective of the vector that is returned.
The practical failure was worse than a suboptimal start. The single optimized
path ran into a degenerate corner in which the numerical integral of the
density underflows, moment_theoretical() falls back to its closed-form
Kumaraswamy moment, and the optimizer is rewarded for parameters whose real
moments are nothing like the sample's. In 15 of 168 sweep cases the returned
vector scored the maximum possible objective of exactly 3.0 -- all five
relative moment errors equal to 1:
set.seed(3050); x <- rgkw(50, 2, 3, 1.5, 2, 0.8) sample mean 0.296299
before alpha = 0.100000 (pinned at the lower bound), beta = 9.874503,
gamma = 0.574504, delta = 0.131984, lambda = 1.276909
theoretical mean 1.401418e-06, objective 3.000000
after alpha = 0.773802, beta = 3.322900, gamma = 3.662037,
delta = 0.670337, lambda = 1.437505, objective 4.830e-07
Downstream the damage reached the fits themselves. Starting optim() from that
corner, 8 of 10 GKw samples of size 500 converged to a positive negative
log-likelihood -- around +4100 to +4900 where the correct region is near -300:
nll from old start nll from new start
gkw seed 1 4117.357 -295.479
gkw seed 4 4829.810 -297.731
gkw seed 10 4855.630 -284.388
Verified over 168 cases (seven families x n in {50, 200, 1000} x 8 seeds) at
the default n_starts = 5: 77 improved, 91 unchanged, none worse, worst ratio
exactly 1.000000. The 15 cases at the maximum objective of 3.0 fell to none,
and the largest objective over the sweep fell from 3.0 to 5.876e-04. Over 70
MLE fits driven from the two sets of starting values, the largest improvement
in the attained negative log-likelihood was 5230.4 nats and the largest
regression 0.09 nats, on a kkw sample that settled in a neighbouring local
optimum; mean relative parameter error fell from 0.507 to 0.409. The full test
suite is unchanged at 0 failures.
Cost: n_starts now buys what it claims, so it also costs what it claims. At
the default n_starts = 5 a GKw call goes from 0.058s to 0.27s (n = 300); at
n_starts = 1000, from 0.09s to 50s. The default is unchanged. Four fixed,
family-specific starting points are always used, so n_starts below 4 still
behaves as 4; this is now documented rather than silently true.
gkwgetstartvalues() truncated out-of-support data without saying so
(gkwinit.cpp): every observation was clamped into [1e-10, 1 - 1e-10] in
silence. Truncation moves every sample moment and therefore every estimate the
function returns, and its commonest cause -- data on a percentage or 0-100
scale -- is exactly the case a caller needs to be told about. On that input the
function did not fail; it answered, and the answer was both wrong and
unremarkable-looking:
set.seed(1); y <- rkw(300, 2, 3)
gkwgetstartvalues(y, "kw") alpha 2.137549 beta 3.506511
gkwgetstartvalues(c(y, 5, -3), "kw") alpha 2.047567 beta 3.251306
gkwgetstartvalues(y * 100, "kw") alpha 50.000000 beta 50.000000
The last line is the whole problem in one row: a sample handed over on a 0-100
scale came back with both parameters pinned at the upper edge of the parameter
box, with nothing to distinguish it from a fit. tryCatch(..., warning = )
caught no condition in any of the three calls.
Observations outside the open interval (0,1) -- exact 0 and exact 1 included,
matching the support ll*() enforces since the fix earlier in this release --
now raise a warning naming how many were truncated and the range they spanned:
gkwgetstartvalues: 300 of 300 observations lie outside the open interval (0,1)
(observed range [6.6169, 89.7703]) and were clamped to it; the estimates below
are those of the clamped sample. Data on a percentage or 0-100 scale must be
rescaled before use.
The clamp itself is kept, so a single boundary observation still does not abort
a fit and no existing call changes its return value; only the silence is
removed. NA and non-finite values continue to be dropped without a warning,
which the @param x entry now states.
Verified: the warning fires on c(y, 5, -3), on a lone exact 0, on a lone
exact 1 and on y * 100, reporting 2, 1, 1 and 300 offenders respectively with
the correct observed range in each; it does not fire on the clean sample, nor
on samples carrying NA or Inf. Every returned vector is unchanged. The full
test suite, whose data are all strictly inside the support, still reports 0
warnings.
Bimodality was claimed in four places and the family does not reach it
(R/gkwdist-package.R, README.Rmd, README.md, vignettes/gkwdist.Rmd):
the overview listed bimodality among the shapes the GKw accommodates, the
"Advantages" list said "bimodal, U-shaped, bathtub", and both shape-selection
tables offered GKw for "Bimodal or U-shaped" data.
No density with two interior modes was found in roughly 330,000 parameter vectors:
search draws 2+ interior modes
structured grid, 7 x 7 x 6 x 5 x 6 8,820 0
log-uniform (0.05, 30)^5, 4001-point grid 200,000 0
log-uniform (1e-3, 300)^5, 3001-point grid 120,000 0
The third sweep first flagged 644 candidates, every one of them with \eqn{\alpha > 160}. At that size \eqn{x^{\alpha}} underflows for all but a sliver next to 1, so a grid uniform in \eqn{x} cannot resolve the density. The extra peaks are jitter: for \eqn{(\alpha, \beta, \gamma, \delta, \lambda) =} (225.7, 70.1, 0.041, 8.04, 2.48) their count grows with resolution -- 2, then 5, then 55 as the grid goes 3001, 30001, 300001 points -- at heights some thirty orders of magnitude below the mode at 19.23. On a grid uniform in \eqn{x^{\alpha}} each flagged case has at most one interior mode.
Every shape observed was monotone, unimodal or U-shaped. A second peak appears
only as a divergence at a boundary, governed by the exponents
\eqn{\alpha\gamma\lambda - 1} and \eqn{\beta(\delta+1) - 1} already
documented in dgkw(). U-shapes and bathtubs are real and those claims stand;
only bimodality is removed. The overview now states what the family does
instead, so that data with two separated interior modes is sent to a mixture
rather than to a larger member of this family. The claim is reported as what
it is -- a search, not a proof.
The vignette's U-shape recommendation named a condition that does not hold
for EKw (vignettes/gkwdist.Rmd): the row replacing the bimodal one first
read "Kumaraswamy or Exponentiated Kw (\eqn{\alpha, \beta < 1})". For EKw
the left-hand exponent is \eqn{\alpha\lambda - 1}, so
\eqn{\alpha, \beta < 1} is not sufficient: dekw(c(1e-8, 1e-4, 0.1), 0.5, 0.5, 3) is 1.9e-05, 1.9e-03, 8.6e-02 -- rising from 0, not a U. The row now
names Kumaraswamy alone, where \eqn{\alpha, \beta < 1} is exactly right.
Two sub-family constraints in the package overview were mathematically
wrong (R/gkwdist-package.R): the "Distribution Family Hierarchy" block
gave Kumaraswamy as "GKw with \eqn{\gamma = \delta = \lambda = 1}" and the
uniform as "set all shape parameters to 1". In this parameterization
\eqn{\delta} is neutral at 0, not 1 -- it enters through
\eqn{B(\gamma, \delta+1)} and \eqn{z^{\delta}} -- so both statements name a
different distribution than the one they claim:
claim max |difference|
Kw(2.3, 3.1) vs GKw(2.3, 3.1, 1, 1, 1) 0.887
Kw(2.3, 3.1) vs GKw(2.3, 3.1, 1, 0, 1) 0 <- correct
Uniform vs GKw(1, 1, 1, 1, 1) 0.960
Uniform vs GKw(1, 1, 1, 0, 1) 0 <- correct
This is the same off-by-one that the deleted _pkgdown.yml block carried,
and it survived there because the overview was the one place where the
constraints were never checked against the implementation. Every other
statement of a sub-family constraint in the package -- the \itemize list in
dgkw()'s details, the @details of each family, and the desc lines in
_pkgdown.yml -- was verified numerically in this pass and is correct. The
uniform entry now says explicitly why \eqn{\delta = 0} rather than
\eqn{\delta = 1} is the neutral value, since that is what the two wrong
statements had in common.
The @return of all seven densities denied the boundary contract
(R/gkw.R, R/bkw.R, R/kkw.R, R/ekw.R, R/bpmc.R, R/kw.R,
R/beta.R): once the closed boundaries began carrying the limiting density,
the seven blocks still promised 0 -- or -Inf on the log scale -- for x
"outside the interval (0, 1)". No rendered page mentioned the limit at all,
and 0 and 1 lie outside that open interval, so the sentence denied precisely
the two points the change had added:
call documented returned
dgkw(c(0, 1), 0.5, 0.5, 1, 0, 1) 0 0 Inf Inf
dkw(c(0, 1), 2, 1) 0 0 0 2
dbeta_(c(0, 1), 0.5, 0) 0 0 Inf 0.5
dgkw(c(-0.1, 1.1), 2, 3, 1, 0, 1) 0 0 0 0
The last row is the case that really does return 0, and the wording is
narrowed to it: "strictly outside the interval [0, 1]". Each block now states
that x = 0 and x = 1 carry the limiting density, names stats::dbeta as
the base R convention being followed -- with the shape2 = delta + 1 shift
spelled out in dbeta_(), whose parameterization differs -- and records that
the limit is 0, a finite positive value, or Inf according to the
parameters.
Documentation only; no executable code changes and no numerical result is affected.
The out-of-support warning was undocumented in all seven ll*()
(R/gkw.R, R/bkw.R, R/kkw.R, R/ekw.R, R/bpmc.R, R/kw.R,
R/beta.R): the guard that warns 'data' contains values outside (0, 1)
appeared in no help page, so the one signal separating a mis-scaled sample
from a genuine fit failure was invisible to a reader of the documentation.
All seven @return blocks now record it, with the reason: an infinite
objective offers an optimiser no gradient direction to follow.
In the same pass llgkw() loses "returns a large positive value (e.g.,
Inf)", the only one of the seven that did not name Inf exactly. Measured,
all seven return Inf exactly, for an invalid par and for out-of-support
data alike.
Documentation only; no executable code changes and no numerical result is affected.
The Cordeiro & de Castro (2011) citation was truncated in 41 of its 42
appearances (the six family files that cite it -- R/kw.R does not): the
entry existed in two broken
shapes -- 14 ended at the journal name with no volume, pages or full stop,
and 27 ended at a dangling comma. Only dgkw() carried it complete. All 42
now read *Journal of Statistical Computation and Simulation*, *81*(7), 883-898.
In the same pass, Carrasco, Ferrari & Cordeiro (2010) -- the primary source
for the five-parameter distribution this package implements -- is added to
the @references of the seven GKw topics. Their @details credited it in
prose while @references listed only the two secondary sources; the earlier
attribution fix reached the prose and not the list. The six sub-families are
left alone: each defines itself as a special case of GKw and links to
dgkw(), so the attribution reaches them through that link.
References were not linkable outside the package overview (all seven
family files, src/gkwinit.cpp): \doi{} appeared on 7 citations in
R/gkwdist-package.R and nowhere else, so 103 citations across the 50
function topics rendered as plain text. Every citation whose DOI was already
recorded in the package -- Carrasco (2010), Cordeiro & de Castro (2011),
Jones (2009), Kumaraswamy (1980) and McDonald (1984) -- now carries it, taken
verbatim from R/gkwdist-package.R and DESCRIPTION. Every topic that has a
\references section now has at least one DOI.
Nadarajah, Johnson/Kotz/Balakrishnan and Devroye are deliberately left
without one: no DOI for those works is recorded anywhere in this package, and
a DOI is not something to reconstruct from memory. gkwgetstartvalues()'s
Jones citation, the only one still in plain text with no \emph, is brought
into the same form as the other seven, and McDonald's volume number is
italicised to match every other volume number in the package.
The package overview page was generated but unreachable
(R/gkwdist-package.R): _PACKAGE carried @keywords internal, which
removes a topic from the help index and from the pkgdown reference.
docs/reference/gkwdist-package.html was being built and nothing linked to
it, so the whole overview -- family hierarchy, performance notes,
model-selection workflow, four worked estimation examples -- could be reached
only by typing the URL. The keyword is dropped and the topic is added to
_pkgdown.yml under a Package Overview heading, which is also what keeps
pkgdown::check_pkgdown() clean once the topic is no longer internal.
Ten of the thirteen \keyword entries were outside R's controlled
vocabulary (all seven family files): density, cumulative, quantile,
random, likelihood, gradient, hessian, beta, kumaraswamy and
mcdonald are absent from R/doc/KEYWORDS, so they indexed nothing. The
seven that name a role become @family groupings -- density, cumulative
distribution, quantile, random generation, log-likelihood, gradient and
Hessian functions, seven members each -- which roxygen2 renders as
bidirectional "Other ...:" links and pkgdown exposes as concepts. The three
that name a distribution become @concept. distribution and optimize,
which are standard, are kept.
This is the axis the hand-written \seealso blocks did not cover: they link
each function to the siblings of its own distribution, never across
distributions, so nothing led from dgkw() to the other six densities. The
curated blocks are untouched and the generated lists are appended below them.
Only beta, kumaraswamy and mcdonald existed as distribution keywords.
The four families that had none -- generalized kumaraswamy,
beta-kumaraswamy, kumaraswamy-kumaraswamy and exponentiated kumaraswamy -- are given the matching concept, so all 49 topics now carry
one family and one distribution concept.
Every example in the package sat inside \donttest{} (all seven family
files): 51 of 52 topics wrapped their entire @examples block, so
R CMD check without --run-donttest -- the form run locally and in most
CI configurations -- executed no example at all, while --as-cran ran all
9,217 lines regardless. The wrapper bought nothing and hid everything.
Measured on one machine, the 28 d/p/q/r blocks take 0.05 s in total
and the 21 ll/gr/hs blocks take 41 s. The wrapper is removed from the
first group and kept on the second, where the cost is real and where the
weakly-identified fits behind the confidence-region entry below make the
examples platform-sensitive. The default check now exercises 1,515 lines of
examples instead of none.
gkwgetstartvalues() had no @seealso (src/gkwinit.cpp): the one
topic whose output is meant to be fed straight into other functions of this
package linked to none of them. It now points at the seven ll*() objectives
it seeds and at stats::optim().
Help page titles followed five competing patterns (all seven family
files): dmc() read "Beta Power Distribution Distribution"; the GKw CDF,
quantile and RNG topics put the role after the distribution ("Generalized
Kumaraswamy Distribution CDF") where the other six put it first; "CDF of the"
and "Cumulative Distribution Function (CDF) of the" both appeared, as did
"Random Generation for" and "Random Number Generation for", and
"Negative Log-Likelihood" took "for the", "for" and "of the" in different
topics. The KKw family was written kkw in all seven of its titles while the
README, _pkgdown.yml and the package overview write KKw.
Sixteen titles are repaired -- the duplicated word, nine against the majority form, and the seven KKw case fixes below -- leaving all 49 on exactly seven patterns, one per role. Titles that abbreviate the distribution in the gradient and Hessian topics are left abbreviated: those are the longest titles in the package and spelling the family out would push them past 120 characters.
The lower-case kkw is corrected in the prose of R/kkw.R as well, 27
occurrences across descriptions, details, @return blocks, @seealso
labels and example plot titles. Three kinds of kkw are deliberately left
alone, because they are identifiers rather than the name of the
distribution: the file references src/kkw.cpp, the exported function names
(dkkw, pkkw, qkkw, rkkw, llkkw, grkkw, hskkw), and the
family = "kkw" argument value that users pass to
gkwgetstartvalues().
Eight verifications inside @examples were commented out (all seven
family files): six q*() topics computed a round trip
p -> q*() -> p*(), printed both numbers, and left the assertion
abs(p_check - p_recalc) < 1e-9 commented. dgkw() and pgkw() were worse:
each built pdf_beta_check / cdf_beta_check against stats::dbeta() /
stats::pbeta() and then commented out the only line that used it, leaving a
variable computed for nothing. All eight now run and print. Checked before
enabling: the two stats comparisons agree to 4.4e-16 and 3.3e-16, and every
round trip is exact to at worst 1.1e-16.
The intro navbar entry pointed at nothing (_pkgdown.yml,
vignettes/): pkgdown fills "Get started" from vignettes/<package>.Rmd,
and the introductory vignette was into-gkwdist.Rmd -- "into" for "intro" --
so the entry was silently dropped and the published site had no "Get started"
link. The file is renamed gkwdist.Rmd, which fixes the spelling and
activates the entry in one move. Its title and VignetteIndexEntry are
unchanged. Links to articles/into-gkwdist.html and calls to
vignette("into-gkwdist") will no longer resolve.
Half of _pkgdown.yml was a superseded copy of itself (_pkgdown.yml):
131 of its 284 lines were a commented-out earlier reference layout. It was
not a
duplicate -- it carried the nesting relation for each sub-family, which the
active version had dropped -- and two of those relations were wrong, giving
EKw as "GKw with γ = δ = 1" and Kw as "GKw with γ = δ = λ = 1" where δ = 0 is
the neutral value in this parameterization and is what every @details block
and every function default states. The nesting is carried over to the active
desc lines with δ = 0, and the dead block is removed.
No inst/WORDLIST (new file): spelling::spell_check_package() reported
191 words and was therefore unusable as a check. A curated list of 45 -- the
family abbreviations, the cited authors, the institutions, terms such as
digamma, trigamma, unimodality and bimodality, and the three
fragments the <doi:...> markup and the quoted package name leave in
DESCRIPTION -- now leaves the .Rd files and DESCRIPTION reporting
nothing, so a real typo will show up. Nothing consumes the list
automatically: there is no tests/spelling.R and spelling is not in
Suggests, so it serves whoever runs the check by hand.
The author name was rendered two ways (R/gkwdist-package.R): 49 topics
say "Lopes, J. E." and the package overview said "J. E. Lopes". The overview
now matches the rest.
Editorial hedging in the rendered help pages (all seven family files):
29 \seealso entries qualified functions that exist and are exported --
"(if these exist)", "(gradient, if available)", "(other functions for this
parameterization, if they exist)". A CRAN help page should not speculate about
its own package's contents. The qualifiers are removed and the informative
half of each label is kept, so \code{grbkw} (gradient, if available) reads
\code{\link{grbkw}} (gradient).
41 cross-references in those same blocks used \code{} where \link{} was
meant, so they rendered as plain text and the help pages could not be
navigated between. All 50 distinct \link{} targets now resolve to an alias
in the package.
Documentation only; no executable code changes and no numerical result is affected.
gkwgetstartvalues() is deterministic, and now says so (gkwinit.cpp):
the help page described "multiple random starting points" without stating that
the randomness is internal and fixed. The extra starting points come from a
generator seeded with a constant, so set.seed() has no effect on the returned
value and .Random.seed is neither read nor advanced:
set.seed(1); a <- gkwgetstartvalues(x, "gkw", 20)
set.seed(9999); b <- gkwgetstartvalues(x, "gkw", 20)
identical(a, b) TRUE
.Random.seed unchanged across the call TRUE
the caller's next runif(1) unchanged TRUE
This is deliberate and is being documented, not changed. The function is a
method-of-moments estimator whose output seeds optimisers elsewhere in a fit;
two calls on the same data must agree, or every downstream fit would inherit a
dependence on the ambient seed and would silently consume draws the caller did
not ask to spend. The lever for a wider search is n_starts, which since the
fix above is both effective and monotone -- more starts can only lower the
objective -- so a seed argument would add a lottery where a monotone control
already exists. A Determinism paragraph in Details now states all three
facts, the @examples block demonstrates them, and a comment at the generator
in gkwinit.cpp records the intent so the constant seed is not mistaken for an
oversight.
The set.seed(123) opening the example is correct and is kept: it makes the
rbeta() sample on the next line reproducible. Its comment now says which of
the two calls it governs.
Documentation only; no numerical result changes.
Confidence-region examples were undrawable where the observed information
was not positive definite (29 @examples blocks across the seven families):
every one built a confidence region from
eigen(solve(hs*(mle, data))[1:2, 1:2]) and then took
diag(sqrt(eig_decomp$values)), with nothing to guarantee the eigenvalues were
non-negative. solve() of an observed information matrix is a covariance
matrix only where that information is positive definite; optim() reports
convergence = 0 on a flat likelihood ridge without establishing it. When an
eigenvalue came back negative, sqrt() produced NaN, the whole region became
NaN, and plot() aborted with need finite 'xlim' values.
The fit these examples rest on is weakly identified -- the observed information
has a condition number between 4.4e+06 and 1.3e+07 -- so which side of the
boundary it lands on depends on the BLAS and the optimiser's path. The examples
passed on Linux and macOS and failed on Windows under --run-donttest.
eigen() is now called with symmetric = TRUE, which is what a covariance
matrix warrants and which keeps the eigenvalues real and ordered, and the
eigenvalues are clamped at zero, so the region degenerates rather than
vanishing. Reproduced against the indefinite block directly: 500 of 500 ellipse
coordinates were NaN before and none are after, and plot() raises the same
need finite 'xlim' values before and succeeds after.
Documentation only; no executable code in R/ or src/ is changed, and every
numerical result is unaffected.
grmc() gradient formula had inverted digamma signs (R/bpmc.R): the
@details block documented psi(gamma + delta + 1) - psi(gamma) for the gamma
component and psi(gamma + delta + 1) - psi(delta + 1) for delta. Since
d log B(gamma, delta+1) / d gamma = psi(gamma) - psi(gamma + delta + 1), both
signs were reversed, and the documented formula disagreed with the returned
value by two orders of magnitude. R/beta.R documented the opposite sign for
the same quantity. The implementation was correct throughout; only the
documentation is changed. The same block's Hessian entry for
d2l/dgamma ddelta in hsmc() had the same sign reversal.
README example 6 inverted the sign of the observed information matrix
(README.Rmd): hsekw() already returns the Hessian of the negative
log-likelihood, so negating it again produced a negative definite matrix and
printed NaN for every asymptotic standard error. Example 3 of the same README
and the vignettes were already correct.
GKw attribution corrected (R/gkw.R): the main help page credited Cordeiro
& de Castro (2011), which introduces the Kw-G family, rather than Carrasco,
Ferrari & Cordeiro (2010), which introduces the five-parameter generalized
Kumaraswamy distribution implemented here. DESCRIPTION, the README and the
vignettes already cited the latter.
A cheat sheet, generated rather than written (cheatsheet/,
pkgdown/assets/cheatsheet/): two A4 pages covering the family tree, all 49
exported functions, the d/p/q/r contract, the shapes the family
reaches, the maximum-likelihood recipe, the per-family par ordering and the
nested-model map. It is published at /cheatsheet/ on the pkgdown site and
linked from the navbar and the home sidebar.
Nothing on the sheet is typed twice, because a sheet that prints the package
version would be stale at the first release. cheatsheet/build.R fills the
template from the package itself: the version from DESCRIPTION, the logo
from man/figures, the six density curves from the package's own d*
functions, and the function count and matrix from the namespace. A gained or
lost export is a build error rather than a quietly wrong sheet.
Rscript cheatsheet/build.R --check exits non-zero when the committed HTML
differs from a fresh build, and the pkgdown workflow runs it before building
the site, so a release that forgets to regenerate the sheet fails CI instead
of publishing the wrong version number. Both directories are in
.Rbuildignore; neither reaches the tarball.
inst/CITATION named the version and the year by hand and had to be
edited at every release to stay true -- the same staleness the cheat sheet is
built to avoid. Both now come from the package metadata: the version from
meta$Version, the year from the Date/Publication field CRAN adds to the
installed DESCRIPTION, falling back to the current year on a development
install where that field does not exist. citation("gkwdist") follows a
version bump on its own.
metadata citation reads
Version 1.1.6, no Date/Publication "(2026) ... version 1.1.6"
Version 1.0.9, Date/Publication 2025-03-04 "(2025) ... version 1.0.9"
tools:::.check_citation() reports no problems.
The README lost what the vignettes already carry. 777 lines to 221: the
seven mathematical specifications now live only in theory-gkwdist, which
derives them with proofs, and seven of the nine worked examples only in
gkwdist. What stays is the overview, the hierarchy, the function table, two
quick starts, and pointers to the cheat sheet, the vignettes and the
reference index.
RcppArmadillo moved out of Imports. It was listed there only because
R/zzz.R carried @import RcppArmadillo, which put import(RcppArmadillo)
in NAMESPACE. Armadillo is header-only for a client package: everything
gkwdist uses from it is compiled into gkwdist.so at install time through
LinkingTo, where RcppArmadillo already appeared. Loading its R namespace
at run time bought nothing, and R CMD check --as-cran reported
Package in Depends/Imports which should probably only be in LinkingTo: 'RcppArmadillo'. The tag, the NAMESPACE entry and the Imports line are
gone; LinkingTo is untouched, so the build is unchanged.
numDeriv moved from Imports to Suggests. No function in R/ calls
grad() or hessian(); the package's own derivatives are analytic and live
in src/. numDeriv is used only by the test suite, as the independent
reference the analytic gr*() and hs*() routines are checked against, and
by \seealso cross-references in the help pages, which Suggests keeps
valid. Every test that calls it is guarded by skip_if_not_installed(), so
the suite runs to completion without it. Installing gkwdist no longer pulls
numDeriv in.
The utils::globalVariables() registration in R/zzz.R is gone. It
listed 39 names, all 39 of which also appear in the sibling package
gkwreg's own registration -- regression-diagnostic artefacts such as
cook_dist, leverage, linpred and model_label that no distributions
package produces, plus "::", ":::" and "log", which are functions, not
variables. The list had been copied across and never pruned. With the whole
call removed, R CMD check still reports
checking R code for possible problems ... OK, so not one of the 39 was
suppressing a real finding. Removing it restores the check's ability to
notice a genuine undefined global in future.
New .github/workflows/sanitizers.yaml runs the compiled code under
runtime sanitizers. CI had none: no ASan, UBSan, valgrind or rhub job
anywhere. That gap is what let the worst defect of this release through. The
SIGFPE from i % 0 on zero-length input, which killed the R process
outright, produced no diagnostic under -Wall -Wextra -Wformat=2, passed
R CMD check cleanly, and kept the whole five-platform R-CMD-check matrix
green. UBSan names it on the first call, with a stack trace:
gkw.cpp:134:21: runtime error: division by zero
#0 ... in dgkw(...) src/gkw.cpp:134
#1 ... in _gkwdist_dgkw src/RcppExports.cpp:415
The workflow has two jobs: UBSan on stock R with GCC, which links libubsan
into gkwdist.so and needs no instrumented R, and R-hub's clang-asan
container, which adds AddressSanitizer. Both are continue-on-error: true
for now, so a finding reports without blocking a pull request.
The re-exported pipe, %>%, is deprecated and will be removed in a future
release. Nothing changes in this one: it is still exported and still
behaves exactly as before.
gkwdist does not use the pipe anywhere -- not in R/, not in the tests, not
in the vignettes. It re-exports magrittr's operator and nothing else, which
is the sole reason magrittr is a hard dependency, so every installation of
a distributions package pulls in a package it never calls. R has had a native
pipe, |>, since 4.1.0, and users who want %>% can attach it from its own
source.
This is announced a release ahead rather than done now because
export("%>%") is public API: code that reads library(gkwdist) and then
uses %>% without attaching magrittr or a tidyverse package would stop
working the moment the export went away, with could not find function "%>%". If your code depends on gkwdist supplying it, switch now to
library(magrittr) (or any tidyverse package that re-exports it), or to
|>. The change is invisible if you already attach magrittr yourself.
test-return-contract.R -- sweeps every parameter of every family through
d, p, q and r with -1, 0, NA and NaN,
asserting they take the identical route, and asserts separately
that delta = 0 is still accepted -- delta is the one parameter whose
bound is >= 0, and the only place the change could have gone off by one.
Also covers the out-of-support warning in all seven ll*(). Fails 60
assertions and errors on 1 against the preceding commit.
New tests/testthat/test-zero-length-input.R covers all 28 routines with
zero-length data and zero-length parameters, the empty-subset idiom, and the
correspondence with the stats package's convention.
New tests/testthat/test-deep-tail-precision.R pins the subnormal and
large-density regimes against a log-space reference. It fails 37 assertions
against 1.1.5.
New tests/testthat/test-density-near-upper-bound.R pins dgkw() in the band
the old guard rejected, together with the nesting identities and the total
mass. It fails 17 assertions against 1.1.5.
New tests/testthat/test-hessian-degenerate.R pins the degenerate cases of
hsgkw() and checks that healthy parameters keep their finite, symmetric
matrices. It fails 8 assertions against 1.1.5.
New tests/testthat/test-mcdonald-no-clamping.R pins llmc(), grmc() and
hsmc() against a closed-form reference, checks the Beta nesting identity and
the boundedness of the objective, and fits a model end to end. It fails 45
assertions against 1.1.5.
New tests/testthat/test-cdf-log-space.R pins the seven CDFs against
log-space references, checks monotonicity, range, F + S = 1, the nesting
identities and agreement with the integral of the density. It fails 9
assertions against 1.1.5.
New tests/testthat/test-quantile-log-space.R pins the seven quantiles
against closed-form inversions, checks that they stay inside (0,1), that
p(q(u)) recovers u, that the boundary conventions are unchanged and that
the nesting identities hold. It fails 28 assertions against 1.1.5.
New tests/testthat/test-startvalues-contract.R pins the three contracts
gkwgetstartvalues() advertises: that n_starts widens the search and never
worsens the fit, that the estimate reproduces the first sample moment, that
the answer is deterministic and leaves .Random.seed alone, and that
truncating out-of-support data raises a warning naming how many observations
were moved. It fails 25 assertions against 1.1.5 and passes 54 on this branch.
test-mle-performance.R compared each scenario's mean parameter error over
its own converged subset, which penalises the more robust scenario: the
analytical gradient fits datasets the numerical baseline gives up on, and
those are the hard ones. After the generators changed, a gkw run converged
on 4 reps the baseline could not touch, whose mean error was 16.1 against 0.50
on the 23 shared reps, and the reported ratio went from 0.82 to 4.66. The
gradient was in fact the better of the two throughout -- it also reached a
lower negative log-likelihood in 21 of those 23 reps. The comparison is now
paired over the reps where both converged, which is what the data generation
in that file was already written for. The corrected check passes against
1.1.5, against the previous commit and against this one.
New tests/testthat/test-rng-log-space.R checks that every generator stays
inside (0,1), that rbkw() reproduces the closed-form inversion of its own
replayed draw, that set.seed() still reproduces, that the sample passes a
Kolmogorov-Smirnov test against its own CDF, and that simulate-then-fit runs
end to end. It fails 8 assertions against 1.1.5.
New tests/testthat/test-argument-contract.R covers the two surfaces the
suite had never touched: the roughly 218 documented stop() conditions in
the R wrappers, and non-finite input. It asserts every bound on every shape
parameter of every d*/p*/q*/r*, the n guard, the log,
lower.tail and log.p guards, and the par-length and data guards of
every ll*/gr*/hs*; then that NA_real_, NaN, +Inf and -Inf are
accepted without error, return one double per input element and leave their
finite neighbours untouched. 526 assertions, and it fails 11 of them against
1.1.5: qkw(NA) returned 1 there, and q*(NaN) collapsed to NA in six of
the seven families. The log-space quantile inversion in this release fixed
both, and this file pins them.
Six of its tests are marked skip(). They assert the behaviour the
functions should have for non-finite input -- d*(NA) giving NA rather
than 0, p*(+Inf) giving 1 rather than 0, q* outside [0, 1]
giving the NaN its own warning promises, and the log = NA and NA-shape
-parameter holes -- and will fail 86 assertions until that defect is fixed.
They are the specification, written down and executable; the fix removes
the skip() line and nothing else.
tests/testthat/test-derivatives-validation.R had 69 tests where its own
header promises 70: BKw was missing Hessian config 3. Restored.
tests/testthat/test-loglikelihood-functions.R asserted
expect_true(result < 0) on all seven families, commented "Log-likelihood
should be negative". The ll*() functions return the negative log
likelihood, whose sign is not a property of anything -- llkw(c(1,1), .) is
0 on a uniform sample and llkw(c(2,2), .) is +2.1. Each test now asserts
the defining identity, ll*(par, data) == -sum(d*(data, ..., log = TRUE)).
dgkw() returned zero for every input (gkw.cpp, utils.h): the package's
numerical helpers log1mexp() and log1pexp() collide with functions of the same
name in R's public Rmath.h API, which use the opposite convention
(log(1 - exp(-x)) for x >= 0). In translation units where Rmath.h's macro was
active, calls bound to R's version, which returns NaN for the negative arguments
used here; every density evaluation then failed its finiteness guard and returned 0.
The helpers are now named gkw_log1mexp() and gkw_log1pexp(). The sub-family
densities were unaffected, as was llgkw(), which routes through vec_log1mexp().
Log-likelihoods of EKw, KKw and BKw were wrong for data near zero
(ekw.cpp, kkw.cpp, bkw.cpp): these routines clamped v = 1 - x^alpha and
w = 1 - v^beta at 1e-10 instead of working in log space. For small x and
moderate alpha, x^alpha rounds to 1 in double precision and w collapses to
zero, so the clamp replaced log(w) = -53 by log(1e-10) = -23. Deviations
reached 6,100 log-units, which silently corrupts AIC, BIC and likelihood ratio
tests. All three families now use the same gkw_log1mexp() formulation as
gkw.cpp, and their scores and Hessians are expressed as ratios of logarithms.
Mixed second derivatives were zeroed at degenerate parameter values
(bkw.cpp, ekw.cpp, kkw.cpp): guards of the form if (abs(p - 1) > eps)
gated mixed partial derivatives that do not carry the vanishing factor. Because
d2l/dalpha dgamma is obtained by differentiating (gamma-1)*log(w) once in
gamma, the (gamma-1) factor is consumed and the term survives at gamma = 1.
hsbkw() returned 0 where the correct value was 271.12; hsekw() and hskkw()
had the same defect at beta = 1.
grkkw() clamped gradient terms at 1000 (kkw.cpp): arbitrary
std::min(..., 1000.0) caps distorted the score in beta by up to 5%. The same
clamp appeared as effective_delta inside llkkw(), capping the likelihood for
delta > 1000. This is the defect removed from ekw.cpp in 1.1.3, which had
survived here.
grkkw() and hskkw() skipped the z block at delta = 0 (kkw.cpp):
the shortcut omitted sum(log(z)) from dl/ddelta and zeroed
d2l/dalpha ddelta, d2l/dbeta ddelta and d2l/ddelta dlambda, none of which
carry a delta factor. delta = 0 is a valid interior value of the likelihood.
The Beta sub-family rejected delta = 0 (utils.h): check_beta_pars()
required delta > 0, unlike the other five validators. Since the sub-family is
parameterised as Beta(gamma, delta + 1), delta = 0 is the legitimate
Beta(gamma, 1) boundary; dbeta_(), pbeta_(), qbeta_(), rbeta_(),
llbeta(), grbeta() and hsbeta() all returned NA/Inf there.
test-boundary-derivatives.R (new): every gradient component and every
Hessian entry of all seven sub-families is now compared individually against two
independent references, the general GKw routines restricted to the constrained
parameter point and numDeriv Richardson extrapolation, over grids that include
the degenerate values gamma = 1, beta = 1, lambda = 1 and delta = 0 and
samples containing observations near zero.
test-density-correctness.R (new): densities are checked to integrate to one,
to agree with the general GKw density at the constrained parameter point, to match
base R for the Beta and closed-form Kumaraswamy cases, and to be the derivative of
the corresponding distribution function. The previous PDF tests asserted only type,
length, non-negativity and finiteness, all of which a vector of zeros satisfies.
Componentwise maximum relative error over 720 parameter configurations per family,
against the general GKw routines and against numDeriv:
| Family | log-likelihood | gradient | Hessian | |:-------|---------------:|---------:|--------:| | GKw | 1.5e-11 | 1.5e-08 | 1.5e-07 | | BKw | 1.6e-12 | 3.7e-08 | 3.6e-08 | | KKw | 3.1e-13 | 3.5e-08 | 4.9e-08 | | EKw | 5.4e-14 | 7.8e-09 | 2.7e-08 | | Mc | 3.5e-13 | 2.5e-09 | 9.5e-10 | | Kw | 4.1e-15 | 7.2e-10 | 1.2e-09 | | Beta | 7.2e-15 | 4.6e-10 | 8.9e-11 |
The residual gradient and Hessian errors are at the accuracy limit of Richardson extrapolation itself; against the GKw reference all seven families agree to 1e-13.
inst/paper/: the JOSS manuscript was rewritten. It now positions the
package explicitly as the distribution layer of the GKw ecosystem, records that
the split from gkwreg was made at the request of JOSS reviewers during that
package's review, reports the measured validation and timing results in place of
the previous unverified figures, and follows the current JOSS AI disclosure
policy. The bibliography was expanded to 19 entries with DOIs verified against
Crossref.
CONTRIBUTING.md and CODE_OF_CONDUCT.md (new): contribution workflow,
support expectations, governance, and the testing standard numerical
contributions are held to. The contributing guide documents the cross-check that
makes bug reports actionable: every sub-family routine must agree with the
general GKw routine evaluated at the corresponding constrained parameter point.
README: the claim that the C++ routines are "10-50x faster than equivalent
R implementations" was replaced with measured figures. The original benchmark
compared -sum(log(dkw(x, 2, 3))) against llkw(), which is C++ against C++
plus R loop overhead, and gives roughly 3x. The genuine gain is in the
derivatives: the analytical score is about 9x faster than Richardson
extrapolation and the analytical Hessian about 38x faster at n = 20,000.
inst/CITATION: updated to version 1.1.5 and pointed at the CRAN canonical
URL; it had been stale at version 1.0.8.
Test coverage rose from 70.8% to 74.2% of combined R and C++ lines. The largest
single gain is in src/gkw.cpp (42.9% to 71.5%), which reflects that dgkw()
now executes its density computation instead of falling through to its
finiteness guard.
test-mle-performance.R: Added skip_on_cran() to all timing-based benchmark
tests. These tests compare wall-clock times of analytical vs. numerical gradients and
are inherently unreliable on shared/loaded CRAN check machines, causing spurious
ERROR results. The tests remain available for local development.llgkw() invalid parameter return (gkw.cpp): Fixed critical error where the
negative log-likelihood returned R_NegInf (−∞) for invalid parameters instead of
R_PosInf (+∞). Gradient-based MLE optimizers interpret −∞ as a global minimum,
causing them to converge to the invalid boundary rather than the true MLE.
gkwinit.cpp — delta validation (gkwinit.cpp): Fixed internal gkw_pdf()
rejecting delta = 0 (a valid GKw parameter value) due to a strict delta <= 0
check that should have been delta < 0.
gkwinit.cpp — EKw/Kw sub-family PDF mapping (gkwinit.cpp): Fixed
ekw_pdf() and kw_pdf() passing delta = 1 instead of the correct delta = 0
when delegating to gkw_pdf(). EKw and Kw are GKw sub-families with delta = 0,
not delta = 1. This produced wrong starting values for MLE of these families.
hsbkw() — v^(β−1) computation (bkw.cpp): Fixed the Hessian of the BKw
negative log-likelihood returning a wrong value for β < 1. The ternary expression
(beta > 1.0) ? v_beta/v : 1.0 coincidentally produces the correct result for
β = 1 but is wrong for all 0 < β < 1. Replaced with the exact formula
safe_exp((beta - 1.0) * ln_v).
safe_exp() underflow scaling (utils.h): Fixed a systematic 10× error in
the moderate-underflow branch. The previous implementation used
DBL_MIN_SAFE * exp(x − log(DBL_MIN)) where DBL_MIN_SAFE = 10 * DBL_MIN,
yielding 10 * exp(x) instead of exp(x). The fix uses
DBL_MIN * exp(x − log(DBL_MIN)) = exp(x) exactly.
dgkw() silent boundary truncation removed (gkw.cpp): Removed a block that
silently skipped data points within SQRT_EPSILON^(1/α) of 0 or 1, returning
density 0 for those points without warning. The log-space computation handles
near-boundary values correctly without this truncation.
llekw() / grekw() — lambda clamping removed (ekw.cpp): Removed the
arbitrary cap lambda_factor = min(lambda_factor, 1000) applied to gradient and
Hessian terms when λ > 1000. This distorted optimization for large-λ scenarios and
produced incorrect standard errors.
gkwinit.cpp: Removed using namespace Rcpp; at file scope; replaced with
explicit Rcpp:: qualifications. Added NA/NaN filtering before moment computation
to prevent silent corruption when input data contains missing values.
bkw.cpp: Removed spurious try/catch blocks wrapping Rcpp::as<arma::vec>()
conversions in grbkw() and hsbkw(). These conversions cannot throw in this
context and the silent fallback masked type errors.
gkw.cpp / ekw.cpp: Refactored Hessian accumulation to build only the upper
triangle inside the observation loop and symmetrize once afterwards with
arma::symmatu(), eliminating O(n × p²) redundant assignments.
utils.h — vec_safe_pow() UB guard: Added guard preventing undefined
behaviour when casting large y_rounded values (> INT_MAX) to int for
odd-exponent sign detection.
utils.h — vec_safe_pow() SIMD fast path: Added an early-return path
arma::exp(y * arma::log(x)) for the common case (y > 0, all x > 0) that
is fully auto-vectorizable, improving throughput in gradient/Hessian evaluation.
gkw.cpp, bkw.cpp, kkw.cpp, ekw.cpp, kw.cpp, bpmc.cpp, beta_.cpp) by removing old commented-out implementations that were kept for reference..Call() invocations and roxygen examples.Analytical derivatives validation (test-derivatives-validation.R):
gr*) and Hessian (hs*) functionsnumDerivMLE performance benchmarks (test-mle-performance.R):
inst/paper/):
This release represents a comprehensive refactoring of the entire package codebase, focusing on numerical stability, code consistency, and maintainability.
Unified utility functions: Introduced utils.h header providing numerically stable implementations of critical functions:
log1mexp(): Stable computation of log(1 - exp(x)) using Mächler (2012) methodologylog1pexp(): Overflow-protected computation of log(1 + exp(x))safe_log(), safe_exp(), safe_pow(): Protected arithmetic operations with graceful handling of edge casesvec_safe_log, vec_log1mexp, etc.) for efficient array operationsConsistent parameter validation: All distribution families now use dedicated parameter checkers (check_pars(), check_kw_pars(), check_ekw_pars(), etc.) that properly handle NaN, Inf, and boundary conditions.
Complete documentation: All C++ source files now include comprehensive Doxygen-style documentation headers describing:
Fixed critical bug in qgkw(): Corrected logic error where lower_tail transformation was incorrectly applied when log_p = TRUE. The probability is now properly converted to linear scale before tail adjustment.
Fixed gradient calculation in grkkw(): Resolved issue where log_z was not recomputed after clamping z to minimum threshold, causing corrupted gradient values near boundaries.
Fixed Hessian calculation in hsmc(): Corrected sign errors and formula for the lambda component of the Hessian matrix for the Beta-Power/McDonald distribution.
Fixed gradient signs in grmc(): Ensured consistent computation of log-likelihood gradient before negation for optimization.
Eliminated unused variables: Removed declared but unused constants (exp_threshold) and intermediate variables across all distribution files.
Removed redundant calculations: Streamlined computations, notably in pgkw() where logarithm was computed twice for the same quantity.
Simplified parameter recycling: Replaced double-modulo indexing pattern (idx = i % k; vec[idx % vec.n_elem]) with direct single-modulo access (vec[i % vec.n_elem]) in random generation functions.
Standardized function signatures: All distribution functions now follow consistent patterns for parameter order, validation, and return value handling.
.dgkw_cpp, .pgkw_cpp, etc.), providing:
All seven distribution families have been refactored with identical improvements:
| Distribution | Parameters | File |
|--------------|------------|------|
| Generalized Kumaraswamy (GKw) | α, β, γ, δ, λ | gkw.cpp |
| Kumaraswamy-Kumaraswamy (KKw) | α, β, δ, λ | kkw.cpp |
| Beta-Kumaraswamy (BKw) | α, β, γ, δ | bkw.cpp |
| Exponentiated Kumaraswamy (EKw) | α, β, λ | ekw.cpp |
| Beta-Power/McDonald (BP/Mc) | γ, δ, λ | bpmc.cpp |
| Kumaraswamy (Kw) | α, β | kw.cpp |
| Beta (GKw-style) | γ, δ | beta.cpp |
Each family includes: density (d*), distribution (p*), quantile (q*), random generation (r*), negative log-likelihood (ll*), gradient (gr*), and Hessian (hs*) functions.
Special thanks to the thorough code review process that identified subtle numerical issues in edge cases, particularly for extreme parameter values and observations near distribution boundaries.
Enhanced Examples for Likelihood Functions: All ll*, gr*, and hs* functions now include comprehensive examples demonstrating:
Professional Visualization Standards:
Complete Coverage: Enhanced documentation for all distribution families (Kw, EKw, KKw, GKw) covering 2 to 5 parameters
Theoretical References: Documentation cites foundational work by Carrasco et al. (2010), Jones (2009), Kumaraswamy (1980), and standard inference theory from Casella & Berger (2002)
gkwgetstartvalues() Functionfamily parameter to support all distribution families
"gkw", "bkw", "kkw", "ekw", "mc", "kw", "beta"optim() with analytical gradientsll*(), gr*(), and hs*() functions use (par, data) signaturegkwfit() functionAuthors@R field formattingnumDeriv)RcppArmadillo only in LinkingTo)(par, data))\left/\right delimiters in GitHub Markdowntestthat
gkwreg provides regression modeling capabilities