0.1.4 was tagged but never sent to CRAN. Its changes, listed below, reach CRAN with this release.
The check that 0.1.4 put in front of the functions that work a
value at a time missed some arguments: p in Qf,
Tol in VMmin, Plot in mlr,
add in ellipse, GradTol, StepTol,
FnTol and ItnLimit in Optim0, and n
in Rnorm, and so in Rmvn, which hands it to
Rnorm. Given a vector, each still stopped with one of R's
own errors. Every argument that drives an if or a
while is now checked, and n in Rmvn
besides: 142 argument positions across 70 functions, among them
n in trapez1 and N in romb, which
gained an if with the fixes below. The help pages of the
nine functions name the arguments added. Since Plot,
add, Double in Dec2Bin and DecFlag
in LOG are TRUE or FALSE, the error, and their help
pages, now ask for "a single value" rather than "a single
number".
At a size of one, several functions failed. runifLC()
at its default n = 1, SymInv, LDLT and
evJacobi on a 1 x 1 matrix, and romb with
N = 1 stopped with an error; trapez0 on a single
point returned NA; trapez1 with n = 1 returned
three times the trapezoidal value, and NashSVD on a
one-column matrix returned 0 as its singular value. In most, a
2:n or 1:(n-1) that should be empty at n = 1
counted down instead (in trapez1 it added both end points
a second time); NashSVD has no pair of columns to rotate
then, and so never set the singular value. Each now handles the
case, as Grad and Deriv2 already did, and the
romb page no longer asks for N above 2.
At n = 0, Rnorm returned FALSE, and Rgamma,
Rgamma0 and Rbeta returned one deviate for most
shapes, where rnorm(0) gives numeric(0). They all
give numeric(0) now.
SQRT(Inf) never returned: it handed Inf to NORM,
whose halving loop cannot end on Inf, and NORM looped
forever on a negative x as well. SQRT(Inf) is now Inf,
with the warning EXP and LOG give; NORM
stops with a message on a non-finite x and gives a negative
fraction for a negative x.
evJacobi stopped with "object 's' not found" on some
symmetric matrices, diag(c(2, 1)) among them: a test on
the rotation sine ran on a pair that needed no rotation, before
any sine had been set. The test now runs only when a rotation is
made.
VMmin never returned when its search direction right
after resetting the Hessian estimate was not downhill, as at a
starting point where the gradient is zero: a misspelt
SmallXCount kept its loop alive. Optim0 stopped
with "object 'RelStep' not found" whenever it ended on its first
iteration; RelStep and DelF are NA then.
MachEps3 and MachEps4 with MaxIter = 0
stopped with "object 'CurEps' not found"; they return their
first trial value, 0.5 and -0.5, with iteration 0.
Dec2Bin(0) took the logarithm of 0 and returned NaN bits;
it returns all zeros, the pattern of +0, and for -0 the same
with the sign bit set. Run.test returned a p-value of 2 when no
nonzero value was left; it returns NA.
Given NA, EXP, LOG and SQRT stopped with
R's "missing value where TRUE/FALSE needed", and SQRT did
the same for NaN. They now return NA for NA, as exp,
log and sqrt do, and SQRT returns NaN for
NaN with the warning EXP and LOG already give.
mlr no longer calls anyNA, which first appeared in
R 3.1.0, above the R (>= 3.0.0) that DESCRIPTION declares.
Source comments that described methods the code does not use
are corrected. erfccheb, betaiapprox and
gammpapprox forward to erfc, betai and
gammp or gammq; erfc goes through
gammq, not a Cephes rational approximation; and what the
EXP guards keep from overflowing is g*P(y) and
Q(y), not y = g*g.
The Ref line on the betacf, betai, erf,
erfc, gammln, gammp, gammq,
gcf, gser, inverfc, invbetai and
invgammp help pages lacked a full stop and ran into the
sentence after it.
GAMMA and LGAMMA hold the Lanczos coefficients the
GNU Scientific Library uses, which only the gammln page
said; their own pages now say it too. gammpapprox and
betaiapprox are no longer titled as quadrature
approximations, which in this package they never were.
The 0.1.4 entries below are corrected where they were wrong. The
EXP error figures were relative errors in units of
.Machine$double.eps, not ulp, and the shortcut was taken
on every call whose result is finite and nonzero, not on every
call. The old EXP also returned Inf at
log(.Machine$double.xmax) itself. Its NaN began near 2e78
rather than 1e78, came back instead of 0 for very large negative
arguments as well, and was nearly but not quite total past 1e82.
The overflow warning below 0.9 * .Machine$double.xmax is
new in 0.1.4, not kept from before, and the 0.1.4 argument check
did not reach every argument it said it did. Smaller corrections
fix the list of bare titles and, under 0.1.0, which functions
carry an attribution.
EXP was not running its own algorithm. The range
reduction left xexp = x/ln2 unrounded, so the reduced
argument g never grew past about 1e-13, the
abs(g) < sqrt(eps) shortcut was taken on every call whose
result is finite and nonzero, and the rational approximation
under it was dead code there. What came back was
2^(x/ln2) and nothing else, up to 223 eps off in relative
terms (about 444 ulp) at the ends of the range, and Inf at
log(.Machine$double.xmax) itself, where x/ln2
rounds to 1024. xexp is now rounded to the nearest
integer, which is what puts g in the interval
abs(g) <= ln2/2 the coefficients were fitted on. The
shortcut is gone with it: once g can be nonzero, taking
exp(g) as 1 throws away eight digits, and the rational
form is exact at g = 0 anyway. The final scaling is done
in two halves, because 2^xexp on its own overflows to Inf,
and underflows to 0, at exponents where the result is still a
finite nonzero number. Worst relative error over the whole range
is now about 1.5 eps rather than 223, about 2 ulp.
EXP returned NaN instead of Inf for very large arguments,
and NaN instead of 0 for very large negative ones, from about
2e78 in magnitude and for nearly every argument past about 1e82.
The overflow guard sat at 0.9 * .Machine$double.xmax, which
is nowhere near where exp overflows, so such an argument reached a
range reduction that could not reduce it: g grew with x
until the numerator and the denominator of the rational form both
overflowed, and their quotient was NaN. The guard is now at
log(.Machine$double.xmax) going up and at
-1075*ln2 going down, the latter being where exp falls
under half the least subnormal. Overflow now warns and returns
Inf for every argument past the upper limit, where it used to
warn only past 0.9 * .Machine$double.xmax; underflow
returns 0 quietly. Both limits agree with exp to the last
double.
EXP(-Inf) returned -Inf. The branch for it had
been the +Inf branch with the sign flipped, so it handed
the argument back instead of the limit. It returns 0, which is
what exp does.
Given a vector, the functions here that work a value at a time
failed with "the condition has length > 1", the message R
raises for an if whose condition is not a single value.
It names neither the function nor the reason. 130 argument
positions across 64 functions, EXP, LOG,
SQRT and the distribution functions among them, are now
checked up front, and the error names the function and the
argument and points at sapply. Arguments that do accept a vector and
return the right answer for one, such as mu in
Qnorm or x in Dnorm, are untouched. Each
of the 64 help pages now names the arguments its own function
constrains, so the ones that are scalar throughout and the
ones that merely have a scalar control argument, such as
standardize in mlr, can be told apart without
reading the source.
The Chol and LDLT examples lost their last line
in the rendered help page. A percent sign starts a comment in
Rd, so L %*% t(L) was cut back to L. Both are
now escaped and the round trip back to M shows.
Help page titles: invbetai was titled as the inverse of
gammp, NORM, DENORM and ellipRange carried
their bare names,
and seven pairs of pages shared one title (conv and
conv0, deconv and deconv0, gammp and
gammq, trapez0 and trapez1, Run.test
and run.p, Pbinom and Ppois, Qnorm
and Qlnorm). Each title now says what its page is about.
LOG, EXP and SQRT say that current CPUs
provide them at the instruction level, so the R versions are
there to be read. DESCRIPTION says the same of the package as a
whole: it is best studied with its source open.
Help pages brought back in line with the code. erf and
erfc call gammp and gammq, not
erfccheb; gammp and gammq call gser
or gcf; betai calls betacf and
gammln; inverfc goes through the normal quantile;
GetP returns solve(H, g); the attribute set by
MachEps3 and MachEps4 is "iteration"; the
mlr columns are R-Student, COV-Ratio and DFFITS;
NORM returns c(f, n); NashSVD takes any real
matrix; GAMMA returns Inf at zero as well as at the
negative integers; the runifLC recurrence adds k.
Distribution parameters are named for what they are: mean and
standard deviation for the normal and log-normal, location and
scale for t, degrees of freedom for chi-square and F, rate for
the gamma density, distribution and quantile functions, scale
for Rgamma and Rgamma0, mean lambda for Poisson
and rate alpha for Rexp.
Spelling and wording fixed across the help pages, among them Decomposition, Meticulous, Numerical, commutative, array, lengths, chi-square and distribution, together with argument names and function names left over from copied pages.
Every help page now carries a runnable example. Where a base R equivalent exists the example shows both, so the agreement can be seen at a glance.
Rbeta was wrong for every case except Beta(1, 1). It built
the two gamma deviates with the shape and the scale arguments
the other way round, so Rbeta(n, 2, 5) had mean 0.65 where the
distribution has 2/7 = 0.286.
Rgamma divided by bet in its alph == 1 branch where the
other two branches multiply, so that one shape came out with
mean 1/bet instead of bet.
Rgamma0 never scaled the accepted deviate by alph, so every
shape returned the same Exp(1) times bet: mean 1 rather than 3
for Rgamma0(n, 3, 1). It also now requires alph >= 1, which is
the range its rejection step is valid over. Use Rgamma
below that.
simps13 had the Simpson weights 4 and 2 on the wrong
points, which left it converging to the wrong value: 21.099
against a true 21.337 even at n = 60. It is now accurate to
1e-12 there.
Both Simpson rules failed at their smallest legal n, 2 and 3, on a seq() with the wrong sign of step.
The simps38 example asked for n = 10, which is not a
multiple of 3, so the help page showed NULL.
The deconv0 example divided by a spectrum containing a
zero and printed Inf and NaN.
The erfccheb title no longer claims a Chebyshev
approximation, which it stopped using in 0.1.0.
ChkFx left show.error.message holding a list
instead of restoring it. options() returns the previous
settings as a list, and that list was being passed back as the
value. The restore now passes the list whole and runs from
on.exit, so early returns are covered too.
mlr(Plot = TRUE) no longer opens a device of its own
and no longer leaves par() changed. It draws on the
current device and restores the previous settings on exit.
Diagnostics in mlr and deconv go through
message() rather than cat(), so they can be
suppressed and do not write to stdout.
The ellipse examples no longer call dev.new().
Packaging for CRAN: Authors@R replaces the free text
author fields, the obsolete LazyLoad field is gone, and the
reference manual is no longer shipped under inst/doc, since CRAN
builds its own and the copy was 86 percent of the tarball.
Qpois and Qbinom now follow the quantile
convention of base R, the smallest k whose cumulative
probability reaches p. They previously returned the largest k
whose cumulative probability did not exceed p, one below
stats::qpois and stats::qbinom for every p. The
two now agree exactly over a wide parameter grid; see tests/.
Qbinom(p, n, pe) returns n for p = 1, and both
functions compare against a tolerance sized to the accuracy of
Ppois and Pbinom, about 1e-12. Feeding a
cumulative probability straight back in as p can still land on
the other side of a tie, because those two are built from the
incomplete gamma and beta rather than the saddle point
algorithms base R uses.
First release.
Calculus, distribution functions, random variate generation, linear algebra, differential equations and optimization, sized for one semester.
The special functions follow freely licensed reference
implementations, each attributed in its own help page, and the
distribution and quantile functions are built on them.
gammln uses the Lanczos g=7 coefficients of the GNU
Scientific Library. gser, gcf, gammp and
gammq follow the Cephes series and continued fraction as
redistributed in ALGLIB. betacf and betai follow
the continued fraction and symmetry transform of the GNU
Scientific Library. inverfc uses the Cephes normal
quantile ndtri.
erf and erfc come from the incomplete gamma
function through erf(x) = P(1/2, x^2) and erfc(x) = Q(1/2, x^2).
That needs no coefficient table of its own and stays accurate
far into the tail: erfc holds to 1e-13 relative down to
erfc(26) = 5.7e-296, where a rational fit in x alone has long
since lost its significant digits.
invgammp and invbetai are bracketed Newton
iterations on the forward functions, following the strategy of
gsl_cdf_gamma_Pinv and gsl_cdf_beta_Pinv. Above the median the
residual is formed in the upper tail, which avoids the loss of
about eight digits that a lower-tail residual suffers for p
near 1.
gammpapprox, betaiapprox and erfccheb are
thin wrappers over gammp, betai and erfc.
Every distribution function is checked against the base R equivalent over a wide parameter grid. See tests/.