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 0About 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 0The 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 1The 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 -> InfEach 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 -> 2Anything 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 bThe 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 0Only 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.00The 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 0p*(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 -> NaNThe 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 densityThe 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 0Both 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 NaNA 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 = InfEvery 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.626571Both 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-49a 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.9262971512a 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.6310211159R::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.81537306log(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.84Both 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.00Correcting 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.430138079907The 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 = InfValues 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 zerosA 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-218A 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.84The 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-10Every 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-07Downstream 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.388Verified 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.000000The 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 0The third sweep first flagged 644 candidates, every one of them with . At that size underflows for all but a sliver next to 1, so a grid uniform in cannot resolve the density. The extra peaks are jitter: for (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 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 and 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 ()”. For EKw the left-hand
exponent is , so 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 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 ”
and the uniform as “set all shape parameters to 1”. In this
parameterization is neutral at 0, not 1 – it enters
through and – 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 <- correctThis 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 rather than 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 0The 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 TRUEThis 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:415The 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.test-derivatives-validation.R):
gr*) and
Hessian (hs*) functionsnumDerivtest-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