Utilities from 'Seminar fuer Statistik' ETH Zurich
Description:
Useful utilities ['goodies'] from Seminar fuer Statistik
ETH Zurich, some of which were ported from S-plus in the 1990s.
For graphics, have pretty (Log-scale) axes eaxis(), an enhanced
Tukey-Anscombe plot, combining histogram and boxplot,
2d-residual plots, a 'tachoPlot()', pretty arrows, etc. For
robustness, have a robust F test and robust range(). For system
support, notably on Linux, provides 'Sys.*()' functions with
more access to system and CPU information. Finally,
miscellaneous utilities such as simple efficient prime numbers,
integer codes, Duplicated(), toLatex.numeric() and is.whole().
Authors:
Martin Maechler [aut, cre] ,
Werner Stahel [ctb] (Functions: compresid2way(), f.robftest(), last(),
p.scales(), p.dnorm()),
Andreas Ruckstuhl [ctb] (Functions: p.arrows(), p.profileTraces(),
p.res.2x()),
Christian Keller [ctb] (Functions: histBxp(), p.tachoPlot()),
Kjetil Halvorsen [ctb] (Functions: KSd(), ecdf.ksCI()),
Alain Hauser [ctb] (Functions: cairoSwd(), is.whole(),
toLatex.numeric()*),
Christoph Buser [ctb] (to function Duplicated()),
Lorenz Gygax [ctb] (to function p.res.2fact()),
Bill Venables [ctb] (Functions: empty.dimnames(), primes()),
Tony Plate [ctb] (to inv.seq()),
Isabelle Flückiger [ctb],
Marcel Wolbers [ctb],
Markus Keller [ctb],
Sandrine Dudoit [ctb],
Jane Fridlyand [ctb],
Greg Snow [ctb] (to loessDemo()),
Henrik Aa. Nielsen [ctb] (to loessDemo()),
Vincent Carey [ctb],
Ben Bolker [ctb],
Philippe Grosjean [ctb],
Frédéric Ibanez [ctb],
Caterina Savi [ctb],
Charles Geyer [ctb],
Jens Oehlschlägel [ctb]
AsciiToInt returns integer codes in 0:255
for each (one byte) character in strings. ichar is an
alias for it, for old S compatibility.
strcodes implements in R the basic engine for translating
characters to corresponding integer codes.
chars8bit() is the inverse function of
AsciiToint, producing “one byte” characters from integer
codes. Note that it (and hence strcodes() depends on the
locale, see Sys.getlocale().
a vector of (unique) character strings, typically of one
character each.
Details
Only codes in 1:127 make up the ASCII encoding which should be
identical for all R versions, whereas the ‘upper’ half
is often determined from the ISO-8859-1 (aka “ISO-Latin 1)”
encoding, but may well differ, depending on the locale setting, see
also Sys.setlocale.
Note that 0 is no longer allowed since, R does not allow
\0 aka nul characters in a string anymore.
Value
AsciiToInt (and hence ichar) and chars8bit return a
vector of the same length as their argument.
strcodes(x, tab) returns a list of the same
length and names as x with list
components of integer vectors with codes in 1:255.
Author(s)
Martin Maechler, partly in 1991 for S-plus
Examples
chars8bit(65:70)#-> "A" "B" .. "F"
stopifnot(identical(LETTERS, chars8bit(65:90)),
identical(AsciiToInt(LETTERS),65:90))## may only work in ISO-latin1 locale (not in UTF-8):
try( strcodes(c(a="ABC", ch="1234", place ="Zürich")))## in "latin-1" gives {otherwise should give NA instead of 252}:## Not run: $a
[1]656667$ch
[1]49505152$place
[1]9025211410599104## End(Not run)
myloc <- Sys.getlocale()if(.Platform $ OS.type =="unix") withAutoprint({# ''should work'' here
try( Sys.setlocale(locale ="de_CH"))# "try": just in case
strcodes(c(a="ABC", ch="1234", place ="Zürich"))# no NA hopefully
AsciiToInt(chars8bit())# -> 1:255 {if setting latin1 succeeded above}
chars8bit(97:140)
try( Sys.setlocale(locale ="de_CH.utf-8"))# "try": just in case
chars8bit(97:140)## typically looks different than above})## Resetting to original locale .. works "mostly":
lapply(strsplit(strsplit(myloc,";")[[1]],"="),function(cc) try(Sys.setlocale(cc[1], cc[2])))-> .scratch
Sys.getlocale()== myloc # TRUE if we have succeeded to reset it
x <-1e7*(-10:50)
y <- dnorm(x, m=10e7, s=20e7)
plot(x,y)## not really nice, the following is better:## For horizontal y-axis labels, need more space:
op <- par(mar=.1+ c(5,5,4,1))
plot(x,y, axes=FALSE, frame=TRUE)
aX <- axTicks(1); axis(1, at=aX, label= axTexpr(1, aX))## horizontal labels on y-axis:
aY <- axTicks(2); axis(2, at=aY, label= axTexpr(2, aY), las=2)
par(op)### -- only 'x' and using log-scale there:
plot(x,y, xaxt="n", log ="x")
aX <- axTicks(1); axis(1, at=aX, label= axTexpr(1, aX))## Now an "engineer's version" ( more ticks; only label "10 ^ k" ) :
axp <- par("xaxp")#-> powers of 10 *inside* 'usr'
axp[3]<-1# such that only 10^. are labeled
aX <- axTicks(1, axp = axp)
xu <-10^ par("usr")[1:2]
e10 <- c(-1,1)+ round(log10(axp[1:2]))## exponents of 10 *outside* 'usr'
v <- c(outer(1:9, e10[1]:e10[2],function(x,E) x *10^ E))
v <- v[xu[1]<= v & v <= xu[2]]
plot(x,y, xaxt="n", log ="x", main ="engineer's version of x - axis")
axis(1, at = aX, label = axTexpr(1, aX, drop.1=TRUE))# 'default'
axis(1, at = v, label =FALSE, tcl =2/3* par("tcl"))
x <- seq(0,10, by =.1)## for matrix, dataframe, .. first lines include a header line:
capture.and.write( cbind(x, log1p(exp(x))), first =5)## first, *middle* and last :
capture.and.write( cbind(x, x^2, x^3), first =4, middle =3, n.dots=1)
For an analysis of variance or regression with (at least) two factors:
Plot components + residuals for two factors according to Tukey's
“forget-it plot”. Try it!
either an aov object with a formula of the form
y ~ a + b, where a and b are factors,
or such a formula.
data
data frame containing a and b.
fac
the two factors used for plotting. Either column numbers or
names for argument data.
label
logical indicating if levels of factors should be shown
in the plot.
numlabel
logical indicating if effects of factors will be shown
in the plot.
xlab, ylab, main
the usual title components, here
with a non-trivial default constructed from aov and the
component factors used.
col, lty, pch
colors, line types, plotting characters to be used
for plotting [1] positive residuals, [2] negative residuals,
[3] grid, [4] labels. If pch is sufficiently long, it will be used
as the list of individual symbols for plotting the y values.
Details
For a two-way analysis of variance, the plot shows the additive
components of the fits for the two factors by the intersections of a
grid, along with the residuals.
The observed values of the target variable are identical to the
vertical coordinate.
The application of the function has been extended to cover more
complicated models. The components of the fit for two factors are
shown as just described, and the residuals are added. The result is a
“component plus residual” plot for two factors in one display.
Value
Invisibly, a list with components
compy
data.frame containing the component effects of the two
factors, and combined effects plus residual
coef
coefficients: Intercept and effects of the factors
Compute numerical derivatives of f() given observations
(x,y), using cubic smoothing splines with GCV, see
smooth.spline. In other words, estimate f′()
and/or f′′() for the model
numeric vectors of same length, supposedly from a model
y ~ f(x).
xout
abscissa values at which to evaluate the derivatives.
spar.offset
numeric fudge added to the smoothing parameter,
see spl.par below.
deriv
integer in 1:2 indicating which
derivatives are to be computed.
spl.spar
direct smoothing parameter for smooth.spline.
If it is NULL (as per default), the smoothing parameter used
will be spar.offset + sp$spar, where sp$spar is the GCV
estimated smoothing parameter, see smooth.spline.
Details
It is well known that for derivative estimation, the optimal smoothing
parameter is larger (more smoothing) than for the function itself.
spar.offset is really just a fudge offset added to the
smoothing parameter. Note that in R's implementation of
smooth.spline, spar is really on the
logλ scale.
When deriv = 1:2 (as per default), both derivatives are
estimated with the same smoothing parameter which is suboptimal
for the single functions individually. Another possibility is to call
D1D2(*, deriv = k) twice with k = 1 and k = 2 and
use a larger smoothing parameter for the second derivative.
Value
a list with several components,
x
the abscissae values at which the derivative(s) are evaluated.
D1
if deriv contains 1, estimated values of
f′(xi) where xi are the values from xout.
D2
if deriv contains 2, estimated values of f′′(xi).
spar
the smoothing parameter used in the (final)
smooth.spline call.
df
the equivalent degrees of freedom in that
smooth.spline call.
Author(s)
Martin Maechler, in 1992 (for S).
See Also
D2ss which calls smooth.spline twice,
first on y, then on the f′(xi) values;
smooth.spline on which it relies completely.
Examples
set.seed(8840)
x <- runif(100,0,10)
y <- sin(x)+ rnorm(100)/4
op <- par(mfrow = c(2,1))
plot(x,y)
lines(ss <- smooth.spline(x,y), col =4)
str(ss[c("df","spar")])
plot(cos,0,10, ylim = c(-1.5,1.5), lwd=2,
main = expression("Estimating f'() : "~~ frac(d,dx)* sin(x)== cos(x)))
offs <- c(-0.1,0,0.1,0.2,0.3)
i <-1for(off in offs){
d12 <- D1D2(x,y, spar.offset = off)
lines(d12$x, d12$D1, col = i <- i+1)}
legend(2,1.6, c("true cos()",paste("sp.off. = ", format(offs))), lwd=1,
col =1:(1+length(offs)), cex =0.8, bg =NA)
par(op)
Compute the numerical first or 2nd derivatives of f() given
observations (x[i], y ~= f(x[i])).
D1tr is the trivial discrete first derivative
using simple difference ratios, whereas D1ss and D2ss
use cubic smoothing splines (see smooth.spline)
to estimate first or second derivatives, respectively.
D2ss first uses smooth.spline for the first derivative
f′() and then applies the same to the predicted values
f^′(ti) (where ti are the values of
xout) to find f^′′(ti).
Usage
D1tr(y, x =1)
D1ss(x, y, xout = x, spar.offset =0.1384, spl.spar=NULL)
D2ss(x, y, xout = x, spar.offset =0.1384, spl.spar=NULL)
Arguments
x, y
numeric vectors of same length, supposedly from a model
y ~ f(x). For D1tr(), x can have length one
and then gets the meaning of h=Δx.
xout
abscissa values at which to evaluate the derivatives.
spar.offset
numeric fudge added to the smoothing parameter(s),
see spl.par below. Note that the current default is there
for historical reasons only, and we often would recommend to use
spar.offset = 0 instead.
spl.spar
direct smoothing parameter(s) for smooth.spline.
If it is NULL (as per default), the smoothing parameter used
will be spar.offset + sp$spar, where sp$spar is the GCV
estimated smoothing parameter for both smooths, see
smooth.spline.
Details
It is well known that for derivative estimation, the optimal smoothing
parameter is larger (more smoothing needed) than for the function itself.
spar.offset is really just a fudge offset added to the
smoothing parameters. Note that in R's implementation of
smooth.spline, spar is really on the
logλ scale.
Value
D1tr() and D1ss() return a numeric vector of the length
of y or xout, respectively.
D2ss() returns a list with components
x
the abscissae values (= xout) at which the
derivative(s) are evaluated.
y
estimated values of f′′(xi).
spl.spar
numeric vector of length 2, contain the spar
arguments to the two smooth.spline calls.
spar.offset
as specified on input (maybe rep()eated to length 2).
Author(s)
Martin Maechler, in 1992 (for S).
See Also
D1D2 which directly uses the 2nd derivative of
the smoothing spline; smooth.spline.
Examples
## First Derivative --- spar.off = 0 ok "asymptotically" (?)
set.seed(330)
mult.fig(12)for(i in1:12){
x <- runif(500,0,10); y <- sin(x)+ rnorm(500)/4
f1 <- D1ss(x=x,y=y, spar.off=0.0)
plot(x,f1, ylim = range(c(-1,1,f1)))
curve(cos(x), col=3, add=TRUE)}
set.seed(8840)
x <- runif(100,0,10)
y <- sin(x)+ rnorm(100)/4
op <- par(mfrow = c(2,1))
plot(x,y)
lines(ss <- smooth.spline(x,y), col =4)
str(ss[c("df","spar")])
xx <- seq(0,10, len=201)
plot(xx,-sin(xx), type ='l', ylim = c(-1.5,1.5))
title(expression("Estimating f''() : "* frac(d^2,dx^2)* sin(x)==-sin(x)))
offs <- c(0.05,0.1,0.1348,0.2)
i <-1for(off in offs){
d12 <- D2ss(x,y, spar.offset = off)
lines(d12, col = i <- i+1)}
legend(2,1.6, c("true : -sin(x)",paste("sp.off. = ", format(offs))), lwd=1,
col =1:(1+length(offs)), cex =0.8, bg =NA)
par(op)
These functions are provided for compatibility with older versions of
the sfsmisc package only, and may be defunct as soon as of the
next release.
Usage
pmax.sa(scalar, arr)
pmin.sa(scalar, arr)
Arguments
scalar
numeric scalar.
arr
any numeric R object, typically array.
Details
pmax.sa(s, a) and pmin.sa(s, a) return (more-dimensional) arrays.
These have been deprecated, because pmax and
pmin do so too, if the array is used as
first argument.
This function implements a simple Gaussian maximum likelihood
discriminant rule, for diagonal class covariance matrices.
In machine learning lingo, this is called “Naive Bayes” (for
continuous predictors). Note that naive Bayes is more general, as it
models discrete predictors as multinomial, i.e., binary predictor
variables as Binomial / Bernoulli.
Usage
dDA(x, cll, pool =TRUE)## S3 method for class 'dDA'
predict(object, newdata, pool = object$pool,...)## S3 method for class 'dDA'
print(x,...)
diagDA(ls, cll, ts, pool =TRUE)
Arguments
x, ls
learning set data matrix, with rows corresponding to
cases (e.g., mRNA samples) and columns to predictor variables
(e.g., genes).
cll
class labels for learning set, must be consecutive integers.
object
object of class dDA.
ts, newdata
test set (prediction) data matrix, with rows corresponding
to cases and columns to predictor variables.
pool
logical flag. If true (by default), the covariance matrices
are assumed to be constant across classes and the discriminant rule
is linear in the data. Otherwise (pool= FALSE), the
covariance matrices may vary across classes and the discriminant
rule is quadratic in the data.
...
further arguments passed to and from methods.
Value
dDA() returns an object of class dDA for which there are
print and predict methods. The latter
returns the same as diagDA():
diagDA() returns an integer vector of class predictions for the
test set.
Author(s)
Sandrine Dudoit, [email protected] and
Jane Fridlyand, [email protected] originally wrote
stat.diag.da() in CRAN package sma which was modified
for speedup by Martin Maechler [email protected]
who also introduced dDA etc.
References
S. Dudoit, J. Fridlyand, and T. P. Speed. (2000)
Comparison of Discrimination Methods for the Classification of Tumors
Using Gene Expression Data.
(Statistics, UC Berkeley, June 2000, Tech Report #576)
Formally, for every element N=x[i], compute the (vector
of) “digits” A of the baseb
representation of the number N, N=∑k=0MAM−kbk.
Revert such a representation to integers.
Usage
digitsBase(x, base =2, ndigits =1+ floor(1e-9+ log(max(x,1), base)))## S3 method for class 'basedInt'
as.integer(x,...)## S3 method for class 'basedInt'
print(x,...)
as.intBase(x, base =2)
bi2int(xlist, base)
Arguments
x
For digitsBase(): non-negative integer (vector) whose
base base digits are wanted.
For as.intBase(): a list of numeric vectors, a character
vector, or an integer matrix as returned by digitsBase(),
representing digits in base base.
base
integer, at least 2 specifying the base for representation.
ndigits
number of bits/digits to use.
...
potential further arguments passed to methods, notably
print.
xlist
a list of integer vectors with entries
typically in 0:(base-1), such as resulting from
digitsBase().
Value
For digitsBase(), an object, say m, of class
"basedInt" which is basically a (ndigits x n)
matrix where m[,i] corresponds to x[i],
n <- length(x) and attr(m,"base") is the input
base.
as.intBase() and the as.integer method for
basedInt objects return an integer vector.
bi2int() is the low-level workhorse of as.intBase().
Note
Some of these functions existed under names digits and
digits.v in previous versions of the sfsmisc package.
Author(s)
Martin Maechler, Dec 4, 1991 (for S-plus; then called digits.v).
Examples
digitsBase(0:12,8)#-- octal representation
empty.dimnames(digitsBase(0:33,2))# binary## This may be handy for just one number (and default decimal):
digits <-function(n, base =10) as.vector(digitsBase(n, base = base))
digits(128982734)# 1 2 8 9 8 2 7 3 4
digits(128, base =8)# 2 0 0## one way of pretty printing (base <= 10!)
b2ch <-function(db)
noquote(gsub("^0+(.{1,})$"," \\1",
apply(db,2, paste, collapse ="")))
b2ch(digitsBase(0:33,2))#-> 0 1 10 11 100 101 ... 100001
b2ch(digitsBase(0:33,4))#-> 0 1 2 3 10 11 12 13 20 ... 200 201## Hexadecimal:
i <- c(1:20,100:106)
M <- digitsBase(i,16)
hexdig <- c(0:9, LETTERS[1:6])
cM <- hexdig[1+ M]; dim(cM)<- dim(M)
b2ch(cM)#-> 1 2 3 4 5 6 7 8 9 A B C D E F 10 11 ... 6A## IP (Internet Protocol) numbers coding: <n>.<n>.<n>.<n> <--> longinteger
ip_ntoa <-function(n)
apply(digitsBase(n, base =256),2, paste, collapse=".")
ip_ntoa(2130706430+(0:9))# "126.255.255.254" ... "127.0.0.7"## and the inverse:
ip_aton <-function(a)
bi2int(lapply(strsplit(a,".", fixed=TRUE), as.integer),256)
n <-2130706430+(0:9)
head(ip <- ip_ntoa(n))
head(ip_aton(ip))
stopifnot( n == ip_aton(ip_ntoa(n )),
ip == ip_ntoa(ip_aton(ip)))## Inverse of digitsBase() : as.integer method for the "basedInt" class
as.integer(M)## or also as.intBase() working from strings:(cb <- apply(digitsBase(0:33,4),2, paste, collapse =""))##-> "000" "001" ..... "200" "201"
all(0:33== as.intBase(cb, base =4))
Duplicated() generalizes the duplicated method for
vectors, by returning indices of “equivalence classes” for
duplicated entries and returning nomatch (NA by default)
for unique entries.
Note that duplicated() is not TRUE for the first time a
duplicate appears, whereas Duplicated() only marks unique
entries with nomatch (NA).
a vector of values that cannot be compared,
passed to both duplicated() and match().
FALSE is a special value, meaning that all values can be
compared, and may be the only value accepted for methods other than
the default. It will be coerced internally to the same type as x.
fromLast
logical indicating if duplication should be considered
from the reverse side, i.e., the last (or rightmost) of identical
elements would correspond to duplicated=FALSE.
nomatch
passed to match(): the value to be
returned in the case when no match is found. Note that it is
coerced to integer.
Value
an integer vector of the same length as v. Can be used as a
factor, e.g., in split,
tapply, etc.
Author(s)
Christoph Buser and Martin Maechler, Seminar fuer Statistik, ETH
Zurich, Sep.2007
An extended axis() function which labels more
prettily, in particular for log-scale axes.
It makes use of plotmath or (LaTeX) expressions of
the form k×10k for labeling a
log-scaled axis and when otherwise exponential formatting would be
used (see pretty10exp).
numeric vector of (“normalsized”) tick locations; by
default axTicks(side, ..), i.e., the same as
axis() would use.
labels
NULL (default), logical,
character or expression, as in axis();
in addition, if NA, labels = TRUE is passed to
axis(), i.e. pretty10exp is not
used. Use FALSE to suppress any labeling.
log
logical or NULL specifying if log-scale should be
used; the default depends on the current plot's axis.
use.expr
logical specifying if pretty10exp(.) should
be used for constructing labels when they are NULL. The
default is typically good enough, but you may occasionally forceuse.expr = TRUE.
f.smalltcl
factor specifying the lengths of the small ticks in
proportion to the normalsized, labeled ticks.
at.small
locations of small ticks; the default,
NULL, uses small.mult and constructs “smart”
locations.
small.mult
positive integer (or NULL), used when
at.small is NULL to indicate which multiples of
at (typically axTicks()) should be used as
“small ticks”. The default NULL will use 9 in
the log case and a number in 2:5 otherwise.
equidist.at.tol
a small positive number, a tolerance to be used
for checking equidistant at values. Used to be hardwired at
.001 which was seen to be too small; increase it when necessary.
small.args
optional list of further arguments to
the (second) axis() call which draws the small ticks.
draw.between.ticks
(only if log is true): logical indicating
that possible (non-small) ticks between the labeled (via at)
ones should be drawn as well (and possibly also used for at.small
construction), see also between.max.
between.max
(only if log and draw.between.ticks
are true): integer indicating ticks should be drawn (approximately)
between the labeled ones.
outer.at
logical specifying that at.small should also be
constructed outside the at range, but still inside the
corresponding par("usr").
drop.1
logical specifying if 1× should be dropped
from labels, passed to pretty10exp().
sub10
logical, integer (of length 1 or 2) or "10", indicating if
some 10k should be simplified to “traditional”
formats, see pretty10exp.
nintLog
only used in R > 2.13.x, when log is true:
approximate (lower bound on) number of intervals for log scaling.
to be set to axp[3] when axp and at
are not specified, in order to tweak the number of (non-small)
tick marks produced from axTicks(..), notably when
log is true, set n.axp to 1, 2, or 3:
maximal number of at values to be used
effectively. If you don't specify at yourself carefully, it
is recommended to set this to something like 25, but this is
not the default, for back compatibility reasons.
las, ...
arguments passed to (the first) axis
call. Note that the default las = 1 differs from
axis's default las = 0.
lab.type
string, passed to pretty10exp to choose
between default plotmath or LaTeX label format.
lab.sep
separator between mantissa and exponent for LaTeX labels,
see pretty10exp.
x <- lseq(1e-10,0.1, length =201)
plot(x, pt(x, df=3), type ="l", xaxt ="n", log ="x")
eaxis(1)## without small ticks:
eaxis(3, at.small=FALSE, col="blue")## If you like the ticks, but prefer traditional (non-"plotmath") labels:
plot(x, gamma(x), type ="l", log ="x")
eaxis(1, labels=NA)
x <- lseq(.001,0.1, length =1000)
plot(x, sin(1/x)*x, type ="l", xaxt ="n", log ="x")
eaxis(1)
eaxis(3, n.axp =1)# -> xaxp[3] = 1: only 10^j (main) ticks## non- log-scale : draw small ticks, but no "10^k" if not needed:
x <- seq(-100,100, length =1000)
plot(x, sin(x)/x, type ="l", xaxt ="n")
eaxis(1)# default -> {1, 2, 5} * 10^j ticks
eaxis(3, n.axp =2)# -> xaxp[3] := 2 -- approximately two (main) ticks
x <- seq(-1,1, length =1000)
plot(x, sin(x)/x, type ="l", xaxt ="n")
eaxis(1, small.args = list(col="blue"))
x <- x/1000
plot(x,1-sin(x)/x, type ="l", xaxt ="n", yaxt ="n")
eaxis(1)
eaxis(2)## more labels than default:
op <- par(lab=c(10,5,7))
plot(x, sin(x)/x, type ="l", xaxt ="n")
eaxis(1)# maybe (depending on your canvas), there are too many,## in that case, maybe use
plot(x, sin(x)/x, type ="l", xaxt ="n")
eaxis(1, axTicks(1)[c(TRUE,FALSE)])# drop every 2nd label
eaxis(3, labels=FALSE)## ore use 'max.at' which thins as well:
plot(x, sin(x)/x, type ="l", xaxt ="n")
eaxis(1, max.at=6)
par(op)### Answering R-help "How do I show real values on a log10 histogram", 26 Mar 2013## the data:
set.seed(1); summary(x <- rlnorm(100, m =2, sdl =3))## the plot (w/o x-axis) :
r <- hist(log10(x), xaxt ="n", xlab ="x [log scale]")## the nice axis:
axt <- axTicks(1)
eaxis(1, at = axt, labels = pretty10exp(10^axt, drop.1=TRUE))## Additionally demo'ing 'sub10' options:
plot(r, xaxt="n")
eaxis(1, at = axt, labels = pretty10exp(10^axt, drop.1=TRUE, sub10 =2))## or
plot(r, xaxt="n")
eaxis(1, at = axt, labels = pretty10exp(10^axt, drop.1=TRUE, sub10 ="10"))## or
plot(r, xaxt="n")
eaxis(1, at = axt, labels = pretty10exp(10^axt, drop.1=TRUE, sub10 = c(-2,2)))
Plots the empirical (cumulative) distribution function (ECDF) for
univariate data, together with upper and lower simultaneous 95% confidence curves,
computed via Kolmogorov-Smirnov' D, see KSd.
Usage
ecdf.ksCI(x, main =NULL, sub =NULL, xlab = deparse(substitute(x)),
ci.col ="red",...)
Compute the prime factorization(s) of integer(s) n.
Usage
factorize(n, verbose =FALSE)
Arguments
n
vector of integers to factorize.
verbose
logical indicating if some progress information should
be printed.
Details
works via primes, currently in a cheap way, sub-optimal
for large composite n.
Value
A named list of the same length as n,
each element a 2-column matrix with column "p" the prime
factors and column~"m" their respective exponents (or
multiplities), i.e., for a prime number n, the resulting matrix
is cbind(p = n, m = 1).
Construct a “list”, really an environment
typically of functions and optionally other R objects, where the
functions and formulas all
all share the same environment.
Consequently, the functions may call each other.
On technical level, this is just a simple wrapper around
list2env().
ee <- funEnv(f =function(x) g(2*(x+1)),
g =function(y) hh(y+1),
hh =function(u) u^2,
info ="Some Information (not a function)")
ls(ee)# here the same as names(ee)## Check that it works: i.e., that "f sees g" and "g sees hh":
stopifnot(all.equal(ee$f(pi),(2*pi+3)^2))
ee$f(0:4)# [1] 9 25 49 81 121
if(interactive()){## Both calls work :
helppdf(Normal)
helppdf("NegBinomial")}elseif(.Platform$OS.type !="windows"){# batch mode (Windows often too slow for this)
od <- setwd(tempdir())
ff <- helppdf(Normal, viewer=NULL)
stopifnot(file.exists(ff)); print(ff)
setwd(od)# revert to previous dir.}
numeric vector of data for histogram. Missing values
(NAs) are allowed.
nclass
recommendation for the number of classes (i.e., bars) the histogram should
have. The default is a number proportional to the logarithm of the length
of x.
breaks
vector of the break points for the bars of the histogram. The count in the
i-th bar is sum(breaks[i] < x <= breaks[i+1])
except that if include.lowest is TRUE (the default),
the first bar also includes points equal to breaks[1]. If
omitted, evenly-spaced break points are determined from
nclass and the extremes of the data.
probability
logical flag: if TRUE, the histogram will be scaled as a probability
density; the sum of the bar heights times bar widths will equal 1. If
FALSE, the heights of the bars will be counts.
include.lowest
If TRUE (the default), the lowest bar will include data
points equal to the lowest break, otherwise it will act like the
other bars (see the description of the breaks argument).
xlab
character or expression for x axis labeling.
...
additional arguments to barplot. The
hist function uses the function barplot to do
the actual plotting; consequently, arguments to the barplot
function that control shading, etc., can also be given to
hist. See the barplot documentation for arguments
angle, density, col, and inside. Do not
use the space or histo arguments.
width
width of the box relative to the height of the histogram. DEFAULT is
0.2.
boxcol
color of filled box. The default is 3.
medcol
the color of the median line. The special value, NA,
indicates the current plotting color (par("col")). The
default is 2. If boxcol=0 and medcol is not
explicitly specified this is set to the current plotting color
(par("col")).
medlwd
median line width. The special value NA, is used to indicate
the current line width (par("lwd")). The default is 5.
whisklty
whisker line type. The special value NA indicates the
current line type (par("lty")). The default is 2
(dotted line).
staplelty
staple (whisker end cap) line type. The special value NA
indicates the current line type (par("lty")). The default is
1 (solid line).
Graphical parameters (see par) may also
be supplied as arguments to this function.
In addition, the high-level graphics arguments described under
par and the arguments to title may be supplied to this
function.
Details
If include.lowest is FALSE the bottom breakpoint must be
strictly less than the minimum of the data, otherwise (the default) it
must be less than or equal to the minimum of the data. The top
breakpoint must be greater than or equal to the maximum of the data.
This function has been called hist.bxp() for 17 years; in 2012,
the increasingly strong CRAN policies required a new name (which could not
be confused with an S3 method name).
Author(s)
S-Plus: Markus Keller, Christian Keller; port to R in 1990's: Martin Mächler.
lab <-"50 samples from a t distribution with 5 d.f."
mult.fig(2*3, main ="Hist() + Rug() and histBxp(*)")for(i in1:3){
my.sample <- rt(50,5)
hist(my.sample, main=lab); rug(my.sample)# for 50 obs., this is ok, too..
histBxp(my.sample, main=lab)}
x <-1:4
integrate.xy(x, exp(x))
print(exp(4)- exp(1), digits =10)# the true integralfor(n in c(10,20,50,100,200)){
x <- seq(1,4, len = n)
cat(formatC(n,wid=4), formatC(integrate.xy(x, exp(x)), dig =9),"\n")}
This function tests whether a numeric or complex vector
or array consists of whole numbers. The function is.integer
is not appropriate for this since it tests whether the vector is of class
integer (see examples).
integer, numeric, or complex vector or
array to be tested
tolerance
maximal distance to the next whole number
Value
The return value has the same dimension as the argument x: if x
is a vector, the function returns a logical vector of the same length;
if x is a matrix or array, the function returns a logical matrix
or array of the same dimensions. Each entry in the result indicates whether
the corresponding entry in x is whole.
## Create a random array, matrix, vector
set.seed(307)
a <- array(runif(24), dim = c(2,3,4))
a[4:8]<-4:8
m <- matrix(runif(12),3,4)
m[2:4]<-2:4
v <- complex(real = seq(0.5,1.5, by =0.1),
imaginary = seq(2.5,3.5, by =0.1))## Find whole entries
is.whole(a)
is.whole(m)
is.whole(v)## Numbers of class integer are always whole
is.whole(dim(a))
is.whole(length(v))
Generate numeric sequences applying a linear recursion nr.it times.
Usage
iterate.lin.recursion(x, coeff, delta =0, nr.it)
Arguments
x
numeric vector with initial values, i.e., specifying
the beginning of the resulting sequence; must be of length (larger
or) equal to length(coeff).
coeff
coefficient vector of the linear recursion.
delta
numeric scalar added to each term; defaults to 0. If not
zero, determines the linear drift component.
nr.it
integer, number of iterations.
Value
numeric vector, say r, of length n + nr.it, where
n = length(x). Initialized as r[1:n] = x, the recursion
is r[k+1] = sum(coeff * r[(k-m+1):k]), where m = length(coeff).
Note
Depending on the zeroes of the characteristic polynomial of coeff,
there are three cases, of convergence, oszillation and divergence.
KSd(90)
KSd(1:9)# now works
op <- par(mfrow=c(2,1))
plot(KSd,10,150)# nice
abline(v = c(75,85), col ="gray")
plot(KSd,79,81, n =1001)# *very* tiny discontinuity at 80
par(op)
integer indicating how many element are desired. If
positive, return the length.out last elements of x; if
negative, the last length.out elements are dropped.
na.rm
logical indicating if the last non-missing value (if any)
shall be returned. By default (it is FALSE and) the last
elements (whatever its values) are returned.
Value
a vector of length abs(length.out) of last values from x.
Note
This function may eventually be deprecated for the standard R
function tail().
a <-1:4; lett <- letters[1:9]; CH <-"Suisse"
all.equal(list (a, lett),
list_(a, lett))# "names for current but not for target"
str(list(a, lett, CH))# [[1]], [[2]], .. (no names)
str(list_(a, lett, CH))# $a $lett ..
stopifnot(identical( list (a, lett, CH),
unname(L <- list_(a, lett, CH))),
is.list(L), names(L)== c("a","lett","CH"),
identical(lett, L$lett)## etc)## The function is currently defined asfunction(...) `names<-`(list(...), vapply(sys.call()[-1L], as.character,""))
A graphical and interactive demonstration and visualization of how
loess works. By clicking on the graphic, the user
determines the current estimation window which is visualized together
with the weights.
logical; if true, polygon(.., density=..)
will be used to shade off the regions where the weights are zero.
w.symbols
logical indicating if the non-zero weights should be
visualized by circles with radius proportional to inch.sym
and w where w are the weights.
sym.col, w.col, line.col
colors for the symbols, weights and
lines, respectively.
Author(s)
As function loess.demo(), written and posted to S-news, on 27
Sep 2001, by Greg Snow, Brigham Young University,
it was modified by Henrik Aa. Nielsen, IMM, DTU,
and subsequently spiffed up for R by Martin Maechler.
if(dev.interactive()){if(requireNamespace("lattice")){
data("ethanol", package ="lattice")
attach(ethanol)
loessDemo(E,NOx, span=.25)
loessDemo(E,NOx, span=.25, family ="symmetric")
loessDemo(E,NOx, degree=0)# Tricube Kernel estimate}## Artificial Example with one outlier
n2 <-50; x <-1:(1+2*n2)
fx <-(x/10-5)^2
y <- fx +4*rnorm(x)
y[n2+1]<-1e4
loessDemo(x,y, span=1/3, ylim= c(0,1000))# not robust !!
loessDemo(x,y, span=1/3, family ="symm")
loessDemo(x,y, span=1/3, family ="symm", w.symb =FALSE, ylim = c(0,40))
loessDemo(x,y, span=1/3, family ="symm", ylim = c(0,40))## but see warnings() --- there's a "fixup"}
(x <- lseq(1,990, length=21))
plot(x, x^4, type ="b", col =2, log ="xy")if(with(R.version, major >=2&& minor >=1))
plot(x, exp(x), type ="b", col =2, log ="xy")
names the file to which LaTeX commands should be written
envir
a string, the LaTeX environment name; default is
"tabular"; useful maybe "array", or other versions of
tabular environments.
nam.center
character specifying row names should be center;
default "l".
col.center
character (vector) specifying how the columns should
be centered; must have values from c("l","c","r"); defaults
to "c".
append
logical; if FALSE, will destroy the file
file before writing commands to it; otherwise (by default),
simply adds commands at the end of file file.
digits
integer; setting of options(digits=..) for
purpose of number representation.
title
a string, possibly using LaTeX commands, which will span
the columns of the LaTeX matrix
Value
No value is returned. This function, when used correctly,
only writes LaTeX commands to a file.
Author(s)
For S: Vincent Carey [email protected], from a
post on Feb.19, 1991 to S-news. Port to R (and a bit more) by
Martin Maechler [email protected].
See Also
latex in package Hmisc is more flexible
(but may surprise by its auto-printing ..).
Examples
mex <- matrix(c(pi,pi/2,pi/4,exp(1),exp(2),exp(3)),nrow=2, byrow=TRUE,
dimnames = list(c("$\\pi$","$e$"), c("a","b","c")))
mat2tex(mex, file = print(tf <- tempfile("mat",,".tex")),
title="$\\pi, e$, etc.")## The last command produces the file "mat<xyz>.tex" containing##> \begin{tabular} {| l|| c| c| c|}##> \multicolumn{ 4 }{c}{ $\pi, e$, etc. } \\ \hline##> \ & a & b & c \\ \hline \hline##> $\pi$ & 3.14 & 1.57 & 0.785 \\ \hline##> $e$ & 2.72 & 7.39 & 20.1 \\ \hline##> \end{tabular}## Now you have to properly embed the contents of this file## in a LaTeX document -- for example, you will need a## preamble, the \begin{document} statement, etc.## Note that the backslash needs protection in dimnames## or title actions.
mat2tex(mex, stdout(), col.center = c("r","r","c"))
missingCh can be used to test whether a value was specified
as an argument to a function. Very much related to the standard R
function missing, here the argument is given by its
name, a character string.
As missingCh() calls missing(), do consider the
caveats about the latter, see missing.
a (function evaluation) environment, in which
the variable named x is to be “missing”.
Value
a logical indicating if the argument named x is
missing in the function “above”, typically the
caller of missingCh, but see the use of envir in the
vapply example.
tst1 <-function(a, b, dd,...)## does not work an with argument named 'c' !
c(b = missingCh("b"), dd = missingCh("dd"))
tst1(2)#-> both 'b' and 'dd' are missing
tst1(,3,,3)## b dd## FALSE TRUE -- as 'b' is not missing but 'dd' is.
Tst <-function(a,b,cc,dd,EEE,...)
vapply(c("a","b","cc","dd","EEE"), missingCh,NA, envir=environment())
Tst()## TRUE ... TRUE -- as all are missing()
Tst(1,,3)## a b cc dd EEE## FALSE TRUE FALSE TRUE TRUE## ..... .....## as 'a' and 'cc' where not missing()## Formal testing:
stopifnot(tst1(),!tst1(,3,3), Tst(),
Tst(1,,3, b=2, E="bar")== c(0,0,1,0,0))## maybe surprising that this ^^ becomes 'dd' and only 'cc' is missing
Do simple matrix plots, providing an easy interface to
matplot by using a default x variable.
Usage
mpl(mat,...)
p.m(mat,...)
Arguments
mat
numeric matrix.
...
further arguments passed to matplot,
e.g., type, xlab, etc.
Details
p.m(m) use the first column of m as x variable,
whereas mpl(m) uses the integers 1, 2, ..., nrow(m)
as coordinates and rownames(m) as axis labels if possible.
Note
These were really created for playing around with curves etc, and
probably should be deprecated since in concrete examples, using
matplot() directly is more appropriate.
Easy Setup for plotting multiple figures (in a rectangular layout) on
one page. It allows to specify a main title and uses smart
defaults for several par calls.
integer; the number of plot figures you'll want to draw.
mfrow, mfcol
instead of nr.plots: integer(2) vectors
giving the rectangular figure layout for par(mfrow = *),
or par(mfcol=*), respectively. The default is to use
mfrow = n2mfrow(nr.plots).
marP
numeric(4) vector of figure margins to add
(“Plus”) to default mar, see below.
mgp
argument for par(mpg= .) with a smaller
default than usual.
mar
argument for par(mar= .) with a smaller
default than usual, using the marP argument, see above.
oma
argument for par(oma= .), by default for
adding space for the main title if necessary.
main
character. The main title to be used for the whole graphic.
tit.wid
numeric specifying the vertical width to be used for the
main title; note that this is only used for the default value of
oma (s. above).
cex.main
numeric; the character size to be used for the main title.
line.main
numeric; the margin line at which the title is written
(via mtext(main, side=3, outer=TRUE, line = line.main, ....)).
col.main, font.main
color and font for main title, passed to
mtext(), see also par(*).
A utility function which basically calls plot(*, type="n")
and text. To have names or numbers instead of points
in a plot is useful for identifaction, e.g., in a residual plot, see
also TA.plot.
Usage
n.plot(x, y =NULL, nam =NULL, abbr = n >=20|| max(nchar(nam))>=8,
xlab =NULL, ylab =NULL, log ="",
cex = par("cex"), col = par("col"),...)
Arguments
x, y
coordinates at which to plot. If y is missing,
x is used for both, if it's a data.frame,
list, 2-column matrix etc – via xy.coords;
formula do not work.
nam
the labels to plot at each (x,y). Per default, these
taken from the data x and y; case numbers 1:n
are taken if no names are available.
abbr
logical indicating if the nam labels should be
abbreviated – with a sensible default.
xlab, ylab
labels for the x- and y- axis, the latter being empty
by default.
log
character specifying if log scaled axes should be used, see
plot.default.
defines relative positiveness of eigenvalues compared
to largest, default=1e-6.
conv.tol
convergence tolerance for algorithm, default=1.0e-7
posd.tol
tolerance for enforcing positive definiteness, default=1.0e-8
maxits
maximum number of iterations
verbose
logical specifying if convergence monitoring should be
verbose.
Details
This implements the algorithm of Higham (2002), then forces symmetry,
then forces positive definiteness using code from
posdefify. This implementation does not make
use of direct LAPACK access for tuning purposes as in the MATLAB code
of Lucas (2001). The algorithm of Knol DL and ten Berge (1989) (not
implemented here) is more general in (1) that it allows contraints to
fix some rows (and columns) of the matrix and (2) to force the
smallest eigenvalue to have a certain value.
the slightly more flexible nearPD which also
returns a classed matrix (class dpoMatrix).
For new code, nearPD() is really preferred to nearcor(),
which hence is considered deprecated.
Add a horizontal boxplot to the current plot. This is mainly an
auxiliary function for histBxp, since
boxplot(*, horizontal = TRUE, add = TRUE) is usually
much preferable to this.
Displays a series of plots of the profile t function and the likelihood
profile traces for the parameters in a nonlinear regression model that
has been fitted with nls and profiled with
profile.nls.
Usage
p.profileTraces(x, cex =1,
subtitle = paste("t-Profiles and traces of ",
deparse(attr(x,"summary")$formula)))
Arguments
x
an object of class "profile.nls", typically resulting from
profile(nls(.)), see
profile.nls.
I <-8; J <-3; K <-20
xx <- factor(rep(rep(1:I, rep(K,I)),J))
yy <- factor(rep(1:J, rep(I*K,J)))
zz <- rt(I*J*K, df=5)#-- Student t with 5 d.f.
p.res.2fact(xx,yy,zz, restr=4, main="i.i.d. t <- 5 random |.| <= 4")
mtext("p.res.2fact(xx,yy,zz, restr= 4, ..)",
line=1, adj=1, outer=TRUE, cex=1)## Real data
data(warpbreaks)(fm1 <- lm(breaks ~ wool*tension, data = warpbreaks))## call via formula method of p.res.2x():
p.res.2x(~ ., fm1)# is shorter than, but equivalent to## p.res.2x(~ wool + tension, fm1) ## or the direct## with(warpbreaks, p.res.2fact(wool, tension, residuals(fm1)))#### whereas this is "transposed":
p.res.2x(~ tension+wool, fm1)
Plot Residuals, e.g., of a multiple linear regression, against two
(predictor) variables, using positively and negatively oriented line
segments for positive and negative residuals.
This is a (S3) generic function with a default and a
formula method.
Usage
p.res.2x(x,...)## Default S3 method:
p.res.2x(x, y, z, restricted, size =1, slwd =1, scol =2:3,
xlab =NULL, ylab =NULL, main =NULL,
xlim = range(x), ylim = range(y),...)## S3 method for class 'formula'
p.res.2x(x =~., data, main = deparse(substitute(data)),
xlab =NULL, ylab =NULL,...)
Arguments
x, y
numeric vectors of the same length specifying 2
covariates. For the formula method, x is a formula.
z
numeric vector of same length as x and y,
typically residuals.
restricted
positive value which truncates the size. The
corresponding symbols are marked by stars.
size
the symbols are scaled so that size is the size of
the largest symbol in cm.
slwd, scol
line width and color(s) for the residual
segments. If scol has length 2 as per default,
the two colors are used for positive and negative z values,
respectively.
xlab, ylab, main
axis labels, and title see title,
each with a sensible default. To suppress, use, e.g., main = "".
xlim, ylim
the basic x- and y- axis extents, see
plot.default. Note that these will be slightly
extended such that segments are not cut off.
...
further arguments passed to plot, or
p.res.2x.default(), respectively.
data
(for the formula method:) a data frame or a fitted
"lm" object.
Details
Each residual zz[i] is visualized as line segment centered at
(xxi,yyi), i=1,…,n, where the
lengths of the segments are proportional to the absolute
values ∥zzi∥.
Positive residuals' line segments have slope +1, and negative ones
slope −1, and scol is used to use different colors for
negative and positive segments.
The formula interface calls p.res.2fact() when
bothx and y are factors.
Author(s)
Andreas Ruckstuhl in June 1991 and Martin Maechler, in 1992, '94, 2003-4.
References
Stahel, W.~A. (2008) Statistische Datenanalyse: Eine
Einführung für Naturwissenschaftler,
5. Auflage, Vieweg, Wiesbaden; Paragraph 13.8.r and 13.8.v.
length of unit (or x and y units) of symbol coordinates in cm.
relsysize
same, as a proportion of the plotting area.
Value
A numeric 2x2 matrix, with rows named x and y, and
columns, named "sy2usr" and "usr2cm" which give the
scale conversion factors from ‘symbol’ (as given) to
‘usr’ coordinates and from these to ‘cm’, respectively.
Author(s)
Werner Stahel, 1990; simplification: M.Maechler, 1993, 2004
See Also
par("usr"), of also ("pin") on which this
is based.
coordinates of points. Numeric vectors of the same length.
Missing values (NAs) are allowed.
angle
numeric vector whose elements give the angles between the
horizontal baseline and the minimum and maximum direction of the
pointer measured clockwise in radians.
size
length of the pointers in cm.
method
string specifying the method to calculate the angle of
the pointer. One of "sensitive", "robust" or
"rank". Only the first two characters are necessary.
The minimum and maximum direction of the pointer corresponds to
min(z) and max(z) if method is "sensitive" or "rank"
and to the upper and lower extreme of z if method is "robust"
(see boxplot or rrange for details). The angle is
proportional to z or rank(z) in case of method="rank".
legend
logical flag: if TRUE (default), a legend giving
the values of the minimum and maximum direction of the pointer is drawn.
show.method
logical flag, defaulting to legend; if true,
the method name is printed.
xlab, ylab
labels for x and y axis; defaults to the
‘expression’ used in the function call.
xlim, ylim
numeric of length 2, the limits for the x and y axis,
respectively; see plot.default.
...
further arguments to plot. Graphical
parameters (see par) may also be supplied as arguments
to this function.
Details
A scatter plot of the variables x and y is plotted. The value of the third
variable z is given by the direction of a pointer (similar to a
tachometer). Observations whose z-coordinate is missing are marked by a dot.
For longer time-series, it is sometimes important to spread the
time-series plots over several subplots.
p.ts(.) does this both automatically, and under manual control.
Actually, this is a generalization of plot.ts
(with different defaults).
a named character vector, with names, the
fields, identical to the names of the list returned
by packageDescription, plus its "file" attribute.
Additionally the resulting vector is of class "Dlist" which
activates a useful print() method.
Note
The file is always returned; not the least that the author
wants to see it quite often as his .libPaths() is
non-trivial and typically longer than 4 entries.
List some system level information about the compiled code library,
typically its dependencies, for R packages with compiled code; for
Unix-alikes or more generally when cmd is installed locally.
character vector of package names of
installedR packages.
cmd
a character string with the name of an OS / system
level program (to be called via system(cmd, ..)) which
gives information about the shared library (of compiled code), also
known as “DLL” (dynamically loadable library) or
“so” ((dynamic) shared object) library. The default,
"ldd" is a standard binary utility on Unix-alike platforms
such as Linux. On macOS, "oTool -L" is used by default.
Details
Note that there seems some language confusion as “DLL” on Windows
is also used for “Dynamic-link Library” and Wikipedia warns about
confusing the two concepts (“dynamically loaded ..” vs
“dynamic-link ..”).
Value
a named list with one entry per package in pkg, the
names being the directory / folder names of the
corresponding pkgs from pkg.
The exact structure of such entries is currently subject to change and
you should not rely on its exact format for now.
# for the example only using standard R packages :
myPkgs <- c("stats","MASS","rpart","Matrix")
pl <- pkgLibs(myPkgs)
pl
stopifnot(exprs ={
is.list(pl)
length(pl)== length(myPkgs)
is.character(pkgD <- names(pl))})## Have seen this failing when a strange development version of "Matrix" was picked up:
try( stopifnot( dir.exists(pkgD)))
For one-dimensional nonparametric regression, plot the data and fitted
values, typically a smooth function, and optionally use segments to
visualize the residuals.
numeric vectors all of the same length, representing
(xi,yi) and fitted (smooth) values y^i.
x will be sorted increasingly if necessary, and yd and
ys accordingly.
Alternatively, ys can be an x-y list (as resulting from
xy.coords) containing fitted values on a
finer grid than the observations x. In that case, the
observational values x[]must be part of the larger
set; seqXtend() may be applied to construct such a set
of abscissa values.
Plots a step function
f(x)= ∑iyi1[ti−1,ti](x), i.e., a piecewise constant function of one variable.
With one argument, plots the empirical cumulative distribution
function.
Evaluate one or several univariate polynomials at several locations,
i.e. compute coef[1] + coef[2]*x + ... + coef[p+1]* x^p
(in the simplest case where x is scalar and coef a vector).
Usage
polyn.eval(coef, x)
Arguments
coef
“numeric” vector or matrix. If a vector, x can be an
array and the result matches x.
If coef is a matrix it specifies several polynomials of the
same degree as rows, x must be a vector, coef[,k] is
for xk−1 and the result
is a matrix of dimension length(x) * nrow(coef).
Note that coef can also be complex or bigrational
(as.bigq(.) from gmp, or arbitrary precision
("mpfr") from Rmpfr, or similar number-like objects
for which basic arithmetic is defined.
x
“numeric” vector or array. Either x or coef must
be a vector.
Details
The stable “Horner rule” is used for evaluation in any case.
When length(coef) == 1L, polyn.eval(coef, x) now returns a
vector of length(x) whereas previously, it just gave the number
coef independent of x.
Value
numeric vector or array, depending on input dimensionalities, see above.
Author(s)
Martin Maechler, ages ago.
See Also
For much more sophisticated handling of polynomials, use the
polynom package, see, e.g., predict.polynomial.
For multivariate polynomials (and also for nice interface to the
orthopolynom package), consider the mpoly package.
a string specifying the method to apply; can be abbreviated.
symmetric
logical, simply passed to eigen (unless
eigen.m is specified); currently, we do not see any reason
for not using TRUE.
eigen.m
the eigen value decomposition of
m, can be specified in case it is already available.
eps.ev
number specifying the tolerance to use, see Details
below.
Details
We form the eigen decomposition
m=VΛV′
where Λ is the
diagonal matrix of eigenvalues, Λj,j=λj, with decreasing eigenvalues λ1≥λ2≥…≥λn.
When the smallest eigenvalue λn are less than
Eps <- eps.ev * abs(lambda[1]), i.e., negative or “almost
zero”, some or all eigenvalues are replaced by positive
(>= Eps) values,
Λ~j,j=λ~j.
Then, m~=VΛ~V′ is computed
and rescaled in order to keep the original diagonal (where that is
>= Eps).
Value
a matrix of the same dimensions and the “same” diagonal
(i.e. diag) as m but with the property to
be positive definite.
Note
As we found out, there are more sophisticated algorithms to solve
this and related problems. See the references and the
nearPD() function in the Matrix package.
We consider nearPD() to also be the successor of this package's nearcor().
Author(s)
Martin Maechler, July 2004
References
Section 4.4.2 of
Gill, P.~E., Murray, W. and Wright, M.~H. (1981)
Practical Optimization, Academic Press.
Cheng, Sheung Hun and Higham, Nick (1998)
A Modified Cholesky Algorithm Based on a Symmetric Indefinite Factorization;
SIAM J. Matrix Anal.\ Appl., 19, 1097–1110.
Knol DL, ten Berge JMF (1989)
Least-squares approximation of an improper correlation matrix by a
proper one.
Psychometrika54, 53–61.
Highham (2002)
Computing the nearest correlation matrix - a problem from finance;
IMA Journal of Numerical Analysis22, 329–343.
Lucas (2001)
Computing nearest covariance and correlation matrices. A thesis
submitted to the University of Manchester for the degree of Master of
Science in the Faculty of Science and Engeneering.
See Also
eigen on which the current methods rely.
nearPD() in the Matrix package.
(Further, the deprecated nearcor() from this package.)
Examples
set.seed(12)
m <- matrix(round(rnorm(25),2),5,5); m <-1+ m + t(m); diag(m)<- diag(m)+4
m
posdefify(m)1000* zapsmall(m - posdefify(m))
Fisher's potato crop data set is of historical interest as an early
example of a multi-factor block design.
Usage
data(potatoes)
Format
A data frame with 64 observations on the following 5 variables.
pos
a factor with levels 1:4.
treat
a factor with 16 levels A to H and
J to Q, i.e., LETTERS[1:17][-9].
nitrogen
a factor specifying the amount of nitrogen
sulfate (NH4), with the four levels 0,1,2,4.
potash
a factor specifying the amount of potassium (K,
‘kalium’) sulfate, with the four levels 0,1,2,4.
yield
a numeric vector giving the yield of potatoes in ...
Source
Bennett, J. H. (1972)
Collected Papers of R. A. Fischer vol.~II, 1925-31;
The University of Adelaide.
References
T.Eden and R. A. Fisher (1929)
Studies in Crop Variation. VI. Experiments on the Response of the
Potato to Potash and Nitrogen.
J. Agricultural Science19, 201–213.
Accessible from Bennett (1972), see above.
Examples
data(potatoes)## See the experimental design:
with(potatoes,{
cat("4 blocks of experiments;","each does every (nitrogen,potash) combination (aka 'treat'ment) once.",'', sep="\n")
print(ftable(table(nitrogen, potash, treat)))
print(ftable(tt <- table(pos,potash,nitrogen)))
tt[cbind(pos,potash,nitrogen)]<- as.character(treat)
cat("The 4 blocks pos = 1, 2, 3, 4:\n")
ftable(tt)})## First plot:
with(potatoes, interaction.plot(potash,nitrogen, response=yield))## ANOVAs:
summary(aov(yield ~ nitrogen * potash + Error(pos), data = potatoes))# "==>" can use simply
summary(aov(yield ~ nitrogen + potash + pos, data = potatoes))# and
summary(aov(yield ~ nitrogen + potash, data = potatoes))
logical indicating if 1× should be
dropped from the resulting expressions.
sub10
logical, "10", a non-negative integer number or
an integer vector of length two, say (k1,k2), indicating if some
10j expressions for j∈J should be formatted
traditionally, notably e.g., 100≡1.
When a (non-negative) number, say k, J={j;j≤k} are all simplified, when a length–2 vector,
J={j;k1≤j≤k2} are.
Special cases: sub10 = TRUE means to use
1 instead of 100 and sub10 = "10" uses both
1 for 100 and 10 for 101; these are short
forms of sub10 = c(0,0) and sub10 = c(0,1) respectively.
digits
number of digits for mantissa (a) construction;
the number of significant digits, see signif.
digits.fuzz
the old deprecated name for digits.
lab.type
a string indicating how the result should look like.
By default, (plotmath-compatible)
expressions are returned. Alternatively,
lab.type = "plotmath" returns LaTeX formatted strings for
labels. (The latter is useful, e.g., when using the tikzDevice
package to generate LaTeX-processed figures.)
lab.sep
character separator between mantissa and exponent for
LaTeX labels; it will be prepended with a backslash,
i.e., ‘"cdot"’ will use ‘"\cdot"’
Value
For the default lab.type = "plotmath",
an expression of the same length as x, typically with elements
of the form a %*% 10 ^ k.
Exceptions are 0 which is kept simple, if drop.1 is
true and a=1, 10 ^ k is used, and if sub10
is not false, a %*% 10 ^ 0 as a, and a %*% 10 ^ k as
as the corresponding formatted number a * 10^k independently of
drop.1.
Otherwise, a character vector of the same length as
x. For lab.type = "latex", currently the only
alternative to the default, these strings are LaTeX (math mode)
compatible strings.
Note
If sub10 is set, it will typically be a small number such as 0,
1, or 2. Setting sub10 = TRUE will be interpreted as
sub10 =1 where resulting exponents k will either be
negative or k≥2.
Author(s)
Martin Maechler; Ben Bolker contributed lab.type = "latex"
and lab.sep.
See Also
axTexpr and eaxis() which build on
pretty10exp(), notably the eaxis() example plots.
The new toLatex.numeric method which gives very similar
results with option scientific = TRUE.
Further, axis, axTicks.
Uses an obvious sieve method (and some care), working with
logical and and integers to be quite fast.
Usage
primes(n, pSeq =NULL)
Arguments
n
a (typically positive integer) number.
pSeq
optionally a vector of primes (2,3,5,...) as if from a
primes() call; must be correct.
The goal is a speedup, but currently we have not found one single
case, where using a non-NULL pSeq is faster.
Details
As the function only uses max(n), n can also be a
vector of numbers.
The famous prime number theorem states that π(n), the
number of primes below n is asymptotically n/log(n) in the sense that limn→∞π(n)⋅log(n)/n∼1.
Equivalently, the inverse of pi(), the n-th prime number
pn is around nlogn; recent results (Pierre Dusart, 1999),
prove that
logn+loglogn−1<npn<logn+loglognforn≥6.
Value
numeric vector of all prime numbers ≤n.
Author(s)
Bill Venables (<= 2001); Martin Maechler gained another 40% speed,
carefully working with logicals and integers.
printTable2() prints a 2-way contingency table “with all
bells and whistles” (currently using German labeling).
margin2table() computes marginals, adds them to the table and
returns a margin2table object the print method for which adds
text decorations (using "-" and "|").
Usage
printTable2(table2, digits =3)
margin2table(x, totName ="sum", name.if.empty=FALSE)## S3 method for class 'margin2table'
print(x, digits =3, quote =FALSE, right =TRUE,...)
Arguments
table2
a matrix with non-negative integer entries, i.e. the
contingency table.
x
a matrix; for print(), the result of margin2table.
digits
Anzahl Dezimalstellen, auf die die Häufigkeiten gerundet
werden sollen.
quote, right
logicals passed to print.default(),
but with different default values.
totName
string to use as row- and column- name if x has
corresponding dimnames.
name.if.empty
logical indicating if the margin “totals”
should be named in any case.
...
further potential arguments, unused currently.
Value
margin2table returns a matrix with added marginals,
i.e., an extra row and column, and is of class "margin2table"
(and "table" still) which has a nice print method.
This is defunct now:
The global DEBUG has been a cheap precursor to R's
options(verbose= .) (or a verbose function argument).
This function prints out its arguments as cat() does,
additionally printing the name of function in which it's been called —
only when a global variable DEBUG exists and is
TRUE.
Closes the PostScript or PDF file
(postscript,pdf), openend by a previous
ps.do (or pdf.latex, or ...) call, using
dev.off, and additionally opens a previewer for that
file, unless the previewer is already up. This almost provides
an ‘interactive’ device (like x11) for
postscript or pdf.
logical, indicating if the postscript or
acrobat reader (e.g., ghostview or acroread or the command
given by command) should be called. By default, find out if
the viewer is already runing on this file and only call it if needed.
command
character, giving a system command for PostScript previewing.
By default, getOption("eps_view") is set to gv -watch -geometry -0+0 -magstep -2 -media BBox -noantialias
which assumes gv (aka ghostview) to be in your OS path.
debug
logical; if TRUE print information during execution.
All functions start a pseudo PostScript or Acrobat preview device, using
postscript or pdf, and further registering
the file name for subsequent calls to pdf.end() or
ps.end().
Usage
pdf.do(file, paper ="default", width =-1, height =-1, onefile =FALSE,
title =NULL, version ="1.4", quiet =FALSE,...)
pdf.latex(file, height =5+ main.space *1.25, width =9.5,
main.space=FALSE, lab.space = main.space,
paper ="special", title =NULL,
lab=c(10,10,7), mgp.lab=c(1.6,0.7,0), mar=c(4,4,0.9,1.1),...)
ps.do(file, width=-1, height=-1, onefile=FALSE, horizontal=FALSE,
title =NULL,...)
ps.latex(file, height =5+ main.space *1.25, width =9.5,
main.space=FALSE, lab.space = main.space,
paper ="special", title =NULL,
lab=c(10,10,7), mgp.lab=c(1.6,0.7,0), mar=c(4,4,0.9,1.1),...)
Arguments
file
character giving the PostScript/PDF file name to be written.
height
device height in inches, height * 2.54 are
cm. The default is 5 plus 1.25 iff main.space.
width
device width in inches; for this and
height, see postscript.
onefile, horizontal
logicals passed to
postscript(..) or pdf(..), most probably
to be left alone.
title
PostScript/PDF (not plot!) title passed to
postscript() or pdf(); by default use a
title with R version and file in it.
version
a string describing the PDF version that will be
required to view the output, see pdf; our (high)
default ensures alpha-transparency.
quiet
logical specifying that some (informative/warning)
messages should not be issued.
main.space
logical; if true, leave space for a main title
(unusual for LaTeX figures!).
lab.space
logical; if true, leave space for x- and y- labels
(by not subtracting from mar).
paper
character (or missing), typically "a4" or
"a4r" in non-America, see postscript. Only
if this is "special" (or missing) are your choices of width
and height completely honored (and this may lead to files that
cannot print on A4) with resizing.
lab
integer of length 3, lab[1:2] are desired number of
tick marks on x- and y- axis, see par(lab=).
mgp.lab
three decreasing numbers determining space for axis
labeling, see par(mgp=), the default is here smaller
than usual.
mar
four numbers, indicating marginal space, see
par(mar=), the default is here smaller than usual.
...
arguments passed to ps.do() or
pdf.do() from
ps.latex / pdf.latex and to ps.options
from ps.do/pdf.do.
Details
ps.latex and pdf.latex have an additional
LaTeX
flavor,
and just differ by some extra par settings from the
*.do siblings: E.g., after ps.do(..)
is called, the graphical parameters c("mar", "mgp", "lab") are
reset (to values that typically are better than the defaults for LaTeX
figures).
Whereas the defaults for paper, width, and heightdiffer between pdf and postscript,
they are set such as to provide very similar functionality, for
the functions ps.do() and pdf.do(); e.g., by default,
both use a full plot on portrait-oriented page of the default paper,
as per getOption("papersize"). pdf.do() sets the default paper to "special"
when both width and height are specified.
Determine the quadrant of planar points, i.e. in which of the four
parts cut by the x- and y- axis the points lie. Zero values
(i.e. points on the axes) are treated as if positive.
Usage
quadrant(x, y=NULL)
Arguments
x, y
numeric vectors of the same length, or x is an x−y
structure and y=NULL, see xy.coords.
Value
numeric vector of same length as x (if that's a vector) with
values in 1:4 indicating the quadrant number of the
corresponding point.
Examples
xy <- as.matrix(expand.grid(x=-7:7, y=-7:7)); rownames(xy)<-NULL(qu <- quadrant(xy))
plot(xy, col = qu+1, main ="quadrant() number", axes =FALSE)
abline(h=0, v=0, col="gray")# the x- and y- axis
text(xy, lab = qu, col = qu+1, adj = c(1.4,0))
logical indicating if the org table has header line (in
the usual "|"-separated org table format).
skip
integer number of initial lines to skip.
encoding
to be used in the main
readLines(file, encoding=encoding) call.
fileEncoding
if file is a file name, i.e., a
character string, and fileEncoding is not the
empty string, file(file, "rt", encoding = fileEncoding) will
be used.
text
instead of file, a character or
string (of a few lines, typically).
quiet
logical to suppress the message
which is signalled when no nrows=* has been specified and the
automatic number of rows is smaller than 95% of the rows / non-header
lines of the file.
...
further arguments passed to read.table.
You should not use encoding (but possibly
fileEncoding!) here, as we do not call
read.table on file (but on a
textConnection).
TODO: It should be easy to extend read.org.table() to also
work for some of the proposed Markdown formats for tables.
Please write to maintainer("sfsmisc") or open a
github issue if you are interested.
t1 <-
"
| a | var2 | C ||---+------+-----||2| may |3.4||7| feb |4.7|
"
d <- read.org.table(text = t1)
d
stopifnot(dim(d)== c(2,3),
identical(names(d), c("a","var2","C")),
d[,"a"]== c(2,7))
Compute the signed relative error componentwise (“vectorized”)
between the target and current vectors,
using the absolute error, i.e., the difference
in case the relative error is not well defined, i.e., when target
is zero or infinite.
relErr():
simply the mean absolute value of the
relative errors between target and current vectors;
typically the “same” as
all.equal.numeric(target, vector, tolerance=0, countEQ=TRUE).
Currently useful only when both vectors are finite.
numeric vector of length() a multiple of
length(target); if an array (incl
matrix), dimensions are preserved; for vectors,
names(target) are preserved.
eps0
non-negative number; values abs(target) < eps0 should
be treated as zero (and hence absolute instead of relative error
be computed). This may be crucial when target is an
"mpfr"-number vector.
Value
relErrV():
a numeric vector of the same length (or array
of the same dimension) as current.
relErr():
a single number.
Author(s)
Martin Maechler,
originally as part of Matrix package's ‘test-tools.R’.
See Also
all.equal.numeric() is similar in spirit but returns TRUE or
string containing the mean relative or absolute error.
Examples
## relErrV() test example: showing how it works fine with {NA, Inf, 0} :
eps <-1e-4*c(-9,-8,-6,-4,0.5,1,5)
target <- c(-1:1,0,0,NA,NaN,Inf,-Inf,Inf,0,Inf,1,-3:3)
current <- c(-1:1,1e-7,NaN,NA,0,Inf,Inf,0,Inf,1,Inf,-3:3+ eps)
cbind(target, current, absE = current-target,
relE = relErrV(target,current))-> M ; M
stopifnot(exprs ={
is.logical(isFr <- is.finite(rF <- M[,"relE"]))
target==current | isFr == is.finite(aF <- M[,"absE"])
identical(aF[!isFr], rF[!isFr])
identical(numeric(), relErrV(numeric(), integer()))# length 0 {used to fail}})
tools::assertError(relErrV(1, numeric()), verbose=TRUE)# no longer allowed## relErr() is pretty simple --- (possibly too simple, currently)
relErr
relErr(target, current)# NA (of course)
all.equal.numeric(target, current)## "'is.NA' value mismatch ..."## comparison after dropping NA's :
hasN <- is.na(target)| is.na(current)
all.equal(target[!hasN], current[!hasN], tolerance=0)# "Mean abs. diff.: Inf"
relErr(target[!hasN], current[!hasN])# NaN (to improve?)## comparison after only keeping cases where both are finite:
finN <- is.finite(target)& is.finite(current)
all.equal(target[finN], current[finN], tol=0)# "Mean abs.d.: 0.000279.."
all.equal(target[finN], current[finN], tol=0, countEQ=TRUE)# " " : 0.000239..
relErr(target[finN], current[finN])# 0.0002392929
Simple constructors of a constant character string from one character,
notably a “blank” string of given string length.
M.M. is now ‘mentally deprecating’ bl.string in
favor of using repChar() in all cases.
With R 3.3.0 (May 2016), the new function
strrep() was introduced; it is faster typically, and more
flexible, e.g. accepting a vector for the 2nd argument.
This (for now informally) deprecates all uses of repChar() and bl.string().
Usage
repChar(char, no)
bl.string(no)
Arguments
char
single character (or arbitrary string).
no
non-negative integer.
Value
One string, i.e., character(1)), for bl.string a
blank string, fulfilling n == nchar(bl.string(n)).
r <- sapply(0:8,function(n) ccat(repChar(" ",n), n))
cbind(r)
repChar("-",4)
repChar("_",6)## it may make sense to a string of more than one character:
repChar("-=- ",6)## show the very simple function definitions:
repChar
bl.string
an integer in {1…26}; the default is
particularly useful.
Details
Note that the default n = 13 makes rotn into
a function that is its own inverse.
Written after having searched for it and found
seqinr::rot13() which was generalized and rendered more
transparently to my eyes.
Value
a character as ch, but with each character (which
belongs to letters or LETTERS
“rotated” by n (positions in the alphabet).
Author(s)
Martin Maechler
See Also
rot2, a completely different rotation (namely in the
plane aka R2).
Examples
rotn(c("ABC","a","b","c"),1)
rotn(c("ABC","a","b","c"),2)
rotn(c("ABC","a","b","c"),26)# rotation by 26 does not change much(ch <- paste("Hello", c("World!","you too")))
rotn(ch)
rotn( rotn(ch ))# rotn(*, 13) is its own inverse
Given a real numbers yi with the particular property that
∑iyi is integer, find integer numbers xi
which are close to yi (∣xi−yi∣<1∀i), and have identical “marginal”
sum, sum(x) == sum(y).
As I found later, the problem is known as “Apportionment Problem”
and it is quite an old problem with several solution methods proposed
historically, but only in 1982, Balinski and Young proved that there
is no method that fulfills three natural desiderata.
Note that the (first) three methods currently available here were all
(re?)-invented by M.Maechler, without any knowledge of the
litterature. At the time of writing, I have not even checked to which
(if any) of the historical methods they match.
character string specifying the algorithm to be used.
Details
Without hindsight, it may be surprising that all three methods give
identical results (in all situations and simulations considered),
notably that the idea of ‘mass shifting’ employed in the
iterative "1greedy" algorithm seems equivalent to the much simpler
idea used in "offset-round".
I am pretty sure that these algorithms solve the Lp
optimization problem, minx∥y−x∥p,
typically for all p∈[1,∞]simultaneously, but have not bothered to find a formal proof.
Value
a numeric vector, say r, of the same length as x, but
with integer values and fulfulling sum(r) == sum(x).
Author(s)
Martin Maechler, November 2007
References
Michel Balinski and H. Peyton Young (1982)
Fair Representation: Meeting the Ideal of One Man, One Vote;
Compute a robust range, i.e. the usual range() as long
as there are no outliers, using the “whisker boundaries” of
boxplot, i.e., boxplot.stats.
Usage
rrange(x, range=1, coef =1.5, na.rm =TRUE)
Arguments
x
numeric vector the robust range of which shall be computed.
range
number for S compatibility; 1.5 * range is
equivalent to coef.
coef
numeric multiplication factor definying the outlier
boundary, see ‘Details’ below.
na.rm
logical indicating how NA values should be
handled; they are simply dropped when na.rm = TRUE as by default.
Details
The robust range is really just what boxplot.stats(x,
coef=coef) returns as the whisker boundaries.
This is the most extreme values x[j] still inside median
plus/minus coef * IQR.
Value
numeric vector c(m,M) with m≤M which is (not
strictly) inside range(x) = c(min(x),max(x)).
A more sophisticated robust range for (strongly) asymmetric data can
be derived from the skewness adjusted boxplot statistics
adjboxStats which is a generalization of
boxplot.stats.
Produce a sequence of unique values (sorted increasingly),
containing the initial set of values x.
This can be useful for setting prediction e.g. ranges in nonparametric
regression.
Usage
seqXtend(x, length., method = c("simple","aim","interpolate"),
from =NULL, to =NULL)
Arguments
x
numeric vector.
length.
integer specifying approximately the desired
length() of the result.
method
string specifying the method to be used. The default,
"simple" uses seq(*, length.out = length.) where
"aim" aims a bit better towards the desired final length,
and "interpolate" interpolates evenly inside
each interval [xi,xi+1] in a way to
make all the new intervalls of approximately the same length.
from, to
numbers to be passed to (the default method for)
seq(), defaulting to the minimal and maximal x
value, respectively.
Value
numeric vector of increasing values, of approximate length
length.
(unless length. < length(unique(x)) in which case, the result
is simply sort(unique(x))),
containing the original values of x.
From, r <- seqXtend(x, *), the original values are at
indices ix <- match(x,r), i.e., identical(x, r[ix]).
Note
method = "interpolate" typically gives the best results. Calling
roundfixS, it also need more computational resources
than the other methods.
Author(s)
Martin Maechler
See Also
seq; plotDS can make particularly
good use of seqXtend()
Examples
a <- c(1,2,10,12)
seqXtend(a,12)# --> simply 1:12
seqXtend(a,12,"interp")# ditto
seqXtend(a,12,"aim")# really worse
stopifnot(all.equal(seqXtend(a,12,"interp"),1:12))## for a "general" x, however, "aim" aims better than default
x <- c(1.2,2.4,4.6,9.9)
length(print(seqXtend(x,12)))# 14
length(print(seqXtend(x,12,"aim")))# 12
length(print(seqXtend(x,12,"int")))# 12## "interpolate" is really nice:
xt <- seqXtend(x,100,"interp")
plot(xt, main="seqXtend(*, 100, \"interpol\")")
points(match(x,xt), x, col =2, pch =20)# .... you don't even see that it's not equidistant# whereas the cheap method shows ...
xt2 <- seqXtend(x,100)
plot(xt2, col="blue")
points(match(x,xt2), x, col =2, pch =20)## with "Date" objects
Drng <- as.Date(c("2007-11-10","2012-07-12"))(px <- pretty(Drng, n =16))# say, for the main labels## say, a finer grid, for ticks -- should be almost equidistant
n3 <-3*length(px)
summary(as.numeric(diff(seqXtend(px, n3))))# wildly varying
summary(as.numeric(diff(seqXtend(px, n3,"aim"))))# (ditto)
summary(as.numeric(diff(seqXtend(px, n3,"int"))))# around 30
Collect (and print) information about the current R session and
environment, using sessionInfo() and more mostly
low-level and platform dependent information.
isRshared() is a utility called from sessionInfoX().
NULL (default), TRUE or a
character vector of R package names, whose
packageDescription()s are wanted. No packages by
default, TRUE takes all currently loaded pkgs.
list.libP
a logical indicating if for all
.libPaths entries, the files should be listed via
list.files.
extraR.env
logical indicating if all environment
variables should be recorded which start with "R_" or
"_R_".
x
typically the result of sessionInfoX().
locale
logical, passed to print.sessionInfo()
indicating if the locale information should be printed.
RLIBS
logical indicating if the information about R_LIBS should
be printed.
Renv
logical indicating if the information about R environment
variables should be printed.
shortRversion()## (including the date, typically for an R Core developer)## but this is shorter:(Rver <- shortRversion(date=FALSE))
shortRversion(spaces=FALSE)# e.g. for a file of even directory name
shortRversion(spaces=FALSE, date=FALSE)# even shorter, ditto## If you want even shorter { abbreviate() will remove spaces, too }:
abbreviate(shortRversion(),11)
abbreviate(shortRversion(date=FALSE),13)
sourceAttach(system.file("test-tools-1.R", package="Matrix", mustWork=TRUE))
search()# shows the new "data base" at position 2## look what it contains:
ls.str(pos =2)
potentical further arguments to be passed to
str; str(utils:::str.default) gives useful list.
Value
invisibly (see invisible) a list with
named components matching the pkgs argument. Each of these
components is a named list with one entry per data(.) argument
name. Each entry is a character vector of the names
of all objects, typically only one.
The side effect is, as with str(), to print
everything (via cat) to the console.
str_data("cluster")
str_data("datasets", max=0, give.attr =FALSE)## Filtering (and return value)
dfl <- str_data("datasets", filterFUN=is.data.frame)
str(df.d <- dfl$datasets)## dim() of all those data frames:
t(sapply(unlist(df.d),function(.) dim(get(.))))### Data sets in all attached packages but "datasets" (and stubs):
s <- search()(Apkgs <- sub("^package:",'', s[grep("^package:", s)]))
str_data(Apkgs[!Apkgs %in% c("datasets","stats","base")])
Return information about the Linux hardware, notably
the CPU (the central processor unit) and memory of the
computer R is running on.
This is currently only available for Linux.
These functions exist on other unix-alike platforms, but produce an
error when called.
name of file the lines of which give the CPU info “as
on Linux”
kind
a character string specifying which
kind of memory is desired.
Value
The Sys.*info() functions return a "simple.list",
here basically a named character vector,
(where the names have been filtered through make.names(*,
unique=TRUE) which is of importance for multi-processor or multi-core
CPUs, such that vector can easily be indexed.
Sys.memGB() returns available memory in giga bytes [GB]; Sys.MIPS() returns a number giving an approximation of
the Million Iinstructions Per Second that
the CPU processes (using “bogomips”). This is a performance
measure of the basic non-numeric processing capabilities.
For single-core Linux systems, often about twice the basic clock rate
in “MHz” (as available by Sys.cpuinfo()["cpu.MHz"]); now,
with multicore systems, the result is often around (but smaller than)
2 * #{cores} * clock.rate.
Note
These currently do rely on the Linux ‘/proc/’ file system, and may not
easily be portable to non-Linux environments.
On multi-processor machines, Sys.cpuinfo() contains each field
for each processor (i.e., names(Sys.cpuinfo()) has
duplicated entries).
Conceivably, the bogoMIPS source code is open and available and could
be built into R.
(n.cores <- parallel::detectCores())if(substr(R.version[["os"]],1,5)=="linux"){##-- only on Linux
Sys.cpuinfo()# which is often ugly; this looks much better:
length(Sys.cpu2 <- local({I <- Sys.cpuinfo(); I[!grepl("^flags", names(I))]}))## may still be too much, notably if n.cores > 2:(Sys3 <- Sys.cpu2[!grepl("[.][0-9]+$", names(Sys.cpu2))])
Sys.MIPS()## just the 'bogomips' from above:
Sys.MIPS()/ as.numeric(Sys.cpuinfo()["cpu.MHz"])## ~~ 2 * #{cores} ((no longer))## Available Memory -- can be crucial:
Sys.memGB()#- default "MemTotal"if(Sys.memGB("MemFree")>16)
message("Be happy! You have more than 16 Gigabytes of free memory")}
character strings of "ALL", specifying which
process status fields are desired.
usefile
logical; if true, system writes to a
temporary file and that is scaned subsequently.
ps.cmd
character string, giving the “ps” command name to be used.
verbose
logical ...
warn.multi
logical ...
Details
Use man ps on your respective Unix system, to see what fields are
supported exactly. Unix dialects do differ here, and,
SunOS-Solaris even has more than one ps command...
Value
Note, that Sys.sizes() currently returns two integers which are
“common” to Solaris and Linux.
From a linear (or glm) model fitted, produce the so-called Tukey-Anscombe
plot. Useful (optional) additions include: 0-line, lowess smooth,
2sigma lines, and automatic labeling of observations.
Result of lm(..), aov(..),
glm(..) or a similar object.
fit
fitted values; you probably want the default here.
res
residuals to use. Default: Weighted ("Pearson") residuals if
weights have been used for the model fit.
labels
strings to use as plotting symbols for each point.
Default(NULL): extract observations' names or use its sequence number.
Use, e.g., "*" to get simple * symbols.
main
main title to plot. Default: sophisticated, resulting in
something like "Tukey-Anscombe Plot of : y ~ x" constructed from
lm.res $ call.
xlab
x-axis label for plot.
draw.smooth
logical; if TRUE, draw a lowess smoother
(with automatic smoothing fraction).
show.call
logical; if TRUE, write the "call"ing syntax with
which the fit was done.
show.2sigma
logical; if TRUE, draw horizontal lines at
±2σ where σ is mad(resid).
lo.iter
positive integer, giving the number of lowess
robustness iterations. The default depends on the model and
is 0 for non Gaussian glm's.
lo.cex
character expansion ("cex") for lowess and other
marginal texts.
par0line
a list of arguments (with reasonable defaults) to be passed to
abline(.) when drawing the x-axis, i.e.,
the y=0 line.
parSmooth, parSigma
each a list of arguments (with reasonable
default) for drawing the smooth curve (if draw.smooth is
true), or the horizontal sigma boundaries (if show.2sigma is
true) respectively.
verbose
logical indicating if some construction details should
be reported (print()ed).
...
further graphical parameters are passed to
n.plot(.).
Side Effects
The above mentioned plot is produced on the current graphic device.
Author(s)
Martin Maechler, Seminar fuer Statistik, ETH Zurich, Switzerland;
[email protected]
See Also
plot.lm which also does a QQ normal plot and more.
Examples
data(stackloss)
TA.plot(lm(stack.loss ~ stack.x))
example(airquality)
summary(lmO <- lm(Ozone ~ ., data= airquality))
TA.plot(lmO)
TA.plot(lmO, label ="O")# instead of case numbersif(FALSE){
TA.plot(lm(cost ~ age+type+car.age, claims, weights=number, na.action=na.omit))}##--- for aov(.) : -------------
data(Gun, package ="nlme")
TA.plot( aov(rounds ~ Method + Physique/Team, data = Gun))##--- Not so clear what it means for GLM, but: ------if(require(rpart)){# for the two datasets only
data(solder, package ="rpart")
TA.plot(glm(skips ~ ., data = solder, family = poisson), cex=.6)
data(kyphosis, package ="rpart")
TA.plot(glm(Kyphosis ~ poly(Age,2)+ Start, data=kyphosis, family = binomial),
cex=.75)# smaller title and plotting characters}
For the case of more than two categories or indices (in INDEX),
traditional tapply(*, simplify = TRUE) still returns a
list when an array may seem more useful and natural. This is provided
by tapplySimpl() if the function FUN() is defined such
as to return a vector of the same length in all cases.
Usage
tapplySimpl(X, INDEX, FUN,...)
Arguments
X
an atomic object, typically a vector. All these arguments
are as in tapply() and are passed to tapply(..).
INDEX
list of (typically more than one) factors, each of same
length as X.
FUN
the function to be applied. For the result to be
simplifiable, FUN() must return a vector of always the same
length.
...
optional arguments to FUN.
Value
If the above conditions are satisfied, the list returned from
r <- tapply(X, INDEX, FUN, ...) is simplified into an
array of rank 1+#{indices}, i.e.,
1+length(INDEX); otherwise, tapplySimpl() returns the list
r, i.e., the same as tapply().
## Using tapply() would give a list (with dim() of a matrix);## here we get 3-array:
data(esoph)
with(esoph,{
mima <<- tapplySimpl(ncases/ncontrols, list(agegp, alcgp), range)
stopifnot(dim(mima)== c(2, nlevels(agegp), nlevels(alcgp)))})
aperm(mima)
logical indicating if rug(y) should be
added to each plot. This is too slow for really large sample sizes.
kernels
character vector of kernel names as allowable for the
kernels argument of the standard density function.
from.f, to.f
numeric giving the left and right limit of the
bandwidth scrollbar.
col
color to be used for the density curve.
Details
library(tcltk) must be working, i.e., Tcl/Tk must have been
installed on your platform, and must have been visible during R's
configuration and/or installation.
You can not only choose the bandwidth (the most important parameter),
but also the kernel, and you can zoom in and out (in x-range only).
Value
none.
(How could this be done? tcltk widgets run as separate processes!)
Author(s)
Martin Maechler, building on demo(tkdensity).
Examples
if(dev.interactive(TRUE))## does really not make sense otherwiseif(try(require("tcltk"))){## sometimes (rarely) there, but broken
data(faithful)
tkdensity(faithful $ eruptions)
set.seed(7)if(require("nor1mix"))
tkdensity(rnorMix(1000, MW.nm9), kernels = c("gaussian","epanechnikov"))}
Formats real numbers, possibly in scientific notation, with a given
number of digits after the decimal point. Output can be used in LaTeX
math mode, e.g., for printing numbers in a table, where each number
has to be printed with the same number of digits after the decimal
point, even if the last digits are zeros.
Usage
## S3 method for class 'numeric'
toLatex(object, digits = format.info(object)[2],
scientific = format.info(object)[3]>0, times ="\\cdot",...)
Arguments
object
a numeric vector.
digits
number of digits after the decimal point (for the
mantissa if scientific). The default behaves the same as R's
format().
scientific
logical indicating if scientific notation a *
10^k should be used. The default behaves the same as R's
format().
times
character string indicating the LaTeX symbol to be used for
the ‘times’ sign.
R does not have S' concept of frame = 0, aka ‘session
frame’. These two function were an attempt to provide a portable
way for working with frame 0, particularly when porting code
from S.
logical or integer specifying you want weekday (‘Wochentag’).
0 or FALSE gives no, 1 or TRUE gives a
short and 2 the long version of the day of the week.
Zeit
logical or integer specifying if time ("Zeit") is desired.
0 or FALSE gives no, 1 or TRUE gives a
hours only and 2 hours and minutes.
Value
A string with the current date/time, in the form specified by the arguments.
The C.* are character vector “constants”,
the German ones actually used by u.Datumvonheute.
Author(s)
Caterina Savi, Martin Maechler
See Also
u.date for a similar English version, and
p.datum which plots.
For English month names, etc month.name.
Compute log() only for high values and keep low ones –
antisymmetrically such that u.log(x) is (once) continuously
differentiable, it computes
f(x)=x for ∣x∣≤c and
sign(x)c⋅(1+log(∣x∣/c))
for ∣x∣≥c.
Usage
u.log(x, c =1)
Arguments
x
numeric vector to be transformed.
c
scalar, > 0
Details
Alternatively, the ‘IHS’ (inverse hyperbolic sine)
transform has been proposed, first in more generality as SU()
curves by Johnson(1949); in its simplest form,
f(x)=arsinh(x)=log(x+x2+1),
which is also antisymmetric, continous and once differentiable as our u.log(.).
Value
numeric vector of same length as x.
Author(s)
Martin Maechler, 24 Jan 1995
References
N. L. Johnson (1949).
Systems of Frequency Curves Generated by Methods of Translation,
Biometrika36, pp 149–176.
doi:10.2307/2332539
See Also
Werner Stahel's sophisticated version of John Tukey's “started log” (which was log(x+c)), with concave extension to negative values and adaptive default choice of c;
logst in his CRAN package relevance, or
LogSt in package DescTools.
Examples
curve(u.log,-3,10); abline(h=0, v=0, col ="gray20", lty =3)
curve(1+ log(x),.01, add =TRUE, col="brown")# simple log
curve(u.log(x,2), add =TRUE, col=2)
curve(u.log(x, c=0.4), add =TRUE, col=4)## Compare with IHS = inverse hyperbolic sine == asinh
ihs <-function(x) log(x+sqrt(x^2+1))# == asinh(x) {aka "arsinh(x)" or "sinh^{-1} (x)"}
xI <- c(-Inf,Inf,NA,NaN)
stopifnot(all.equal(xI, asinh(xI)))# but not for ihs():
cbind(xI, asinh=asinh(xI), ihs=ihs(xI))# differs for -Inf
x <- runif(500,0,4); x[100+0:3]<- xI
all.equal(ihs(x), asinh(x))#==> is.NA value mismatch: asinh() is correct {i.e. better!}
curve(u.log,-2,20, n=1000); abline(h=0, v=0, col ="gray20", lty =3)
curve(ihs(x)+1-log(2), add=TRUE, col=adjustcolor(2,1/2), lwd=2)
curve(ihs(x), add=TRUE, col=adjustcolor(4,1/2), lwd=2)## for x >= 0, u.log(x) is nicely between IHS(x) and shifted IHS## a log10-scale version of asinh() {aka "IHS" }: ihs10(x) := asinh(x/2) / ln(10)
ihs10 <-function(x) asinh(x/2)/log(10)
xyaxis <-function() abline(h=0, v=0, col ="gray20", lty =3)
leg3 <-function(x ="right")
legend(x, legend = c(quote(ihs10(x)== asinh(x/2)/log(10)),
quote(log[10](1+x)), quote(log[10](x))),
col=c(1,2,5), bty="n", lwd=2)
curve(asinh(x/2)/log(10),-5,20, n=1000, lwd=2); xyaxis()
curve(log10(1+x), col=2, lwd=2, add=TRUE)
curve(log10( x ), col=5, lwd=2, add=TRUE); leg3()## zoom out and x-log-scale
curve(asinh(x/2)/log(10),.1,100, log="x", n=1000); xyaxis()
curve(log10(1+x), col=2, add=TRUE)
curve(log10( x ), col=5, add=TRUE); leg3("center")
curve(log10(1+x)- ihs10(x),.1,1000, col=2, n=1000, log="x", ylim = c(-1,1)*0.10,
main ="absolute difference", xaxt="n"); xyaxis(); eaxis(1, sub10=1)
curve(log10( x )- ihs10(x), col=4, n=1000, add =TRUE)
curve(abs(1- ihs10(x)/ log10(1+x)),.1,5000, col=2, log ="xy", ylim = c(6e-9,2),
main="|relative error| of approx. ihs10(x) := asinh(x/2)/log(10)", n=1000, axes=FALSE)
eaxis(1, sub10=1); eaxis(2, sub10=TRUE)
curve(abs(1- ihs10(x)/ log10(x)), col=4, n=1000, add =TRUE)
legend("left", legend = c(quote(log[10](1+x)), quote(log[10](x))),
col=c(2,4), bty="n", lwd=1)## Compare with Stahel's version of "started log"## (here, for *vectors* only, and 'base', as Desctools::LogSt();## by MM: "modularized" by providing a threshold-computer function separately:
logst_thrWS <-function(x, mult =1){
lq <- quantile(x[x >0], probs = c(0.25,0.75), na.rm =TRUE, names =FALSE)if(lq[1]== lq[2])
lq[1]<- lq[2]/2
lq[1]^(1+ mult)/lq[2]^mult
}
logst0L <-function(x, calib = x, threshold = thrFUN(calib, mult=mult),
thrFUN = logst_thrWS, mult =1, base =10){## logical index sub-assignment instead of ifelse(): { already in DescTools::LogSt }
res <- x # incl NA's
notNA <-!is.na(sml <-(x <(th <- threshold)))
i1 <- sml & notNA; res[i1]<- log(th, base)+((x[i1]- th)/(th * log(base)))
i2 <-!sml & notNA; res[i2]<- log(x[i2], base)
attr(res,"threshold")<- th
attr(res,"base")<- base
res
}
logst0 <-function(x, calib = x, threshold = thrFUN(calib, mult=mult),
thrFUN = logst_thrWS, mult =1, base =10){## Using pmax.int() instead of logical indexing -- NA's work automatically - even faster
xm <- pmax.int(threshold, x)
res <- log(xm, base)+(x - xm)/(threshold * log(base))
attr(res,"threshold")<- threshold
attr(res,"base")<- base
res
}## u.log() is really using natural log() -- whereas logst() defaults to base=10
curve(u.log,-4,10, n=1000); abline(h=0, v=0, col ="gray20", lty =3); points(-1:1,-1:1, pch=3)
curve(log10(x)+1, add=TRUE, col=adjustcolor("midnightblue",1/2), lwd=4, lty=6)
curve(log10(x), add=TRUE, col=adjustcolor("skyblue3",1/2), lwd=4, lty=7)
curve(logst0(x, threshold=2), add=TRUE, col=adjustcolor("orange",1/2), lwd=2)
curve(logst0(x, threshold=1), add=TRUE, col=adjustcolor(2,1/2), lwd=2)
curve(logst0(x, threshold=1/4), add=TRUE, col=adjustcolor(3,1/2), lwd=2, lty=2)
curve(logst0(x, threshold=1/8), add=TRUE, col=adjustcolor(4,1/2), lwd=2, lty=2)
u.sys() is a convenient wrapper (of system()) to call to
the underlying operating system. The main purpose has been to provide
a function with identical UI both in S-PLUS and R.
MM thinks you shouldn't use this anymore, usually.
Sys.ps.cmd() returns the ‘ps’ (‘process status’)
OS command name (as character string), and is typically
usable on unix alikes only.
Usage
u.sys(..., intern =TRUE)
Sys.ps.cmd()
Arguments
...
any number of strings – which will be
paste()d together and passed to system.
intern
logical – note that the default is reversed from
the one in system().
Author(s)
Martin Maechler
See Also
system, really!;
on non-Windows, Sys.ps() which makes use of Sys.ps.cmd().
Examples
u.sys # shows how simply the function is defined :## Not run: function(..., intern =TRUE)
system(paste(..., sep =""), intern = intern)## End(Not run)# All *running* processes of user [sometimes only R]:
try ( u.sys(Sys.ps.cmd(),"ur"))
Give regularly spaced points on interval [−c,c] with mean 0
(exactly) and variance about 1 (very close for evenn
and larger round.dig). Note that c depends on n.
Usage
unif(n, round.dig =1+ trunc(log10(n)))
Arguments
n
positive integer specifying the number of points desired.
round.dig
integer indicating to how many digits the result is
rounded.
Value
numeric vector of length n, symmetric around 0, hence
with exact mean 0, and variance approximately 1.
Note
It relies on the fact that Var(1,2,...,n)=n(n+1)/12.
logical vector of the same length as x. For the
reversion to work this should select at least all unique values of
x.
need.sort
logical indicating if x is not yet sorted.
Note that this argument exists only for speedup possibility when it
is known, and that it must be set correctly.
Concatenate vector elements or anything using
paste(*, collapse = .).
These are simple short abbreviations I have been using in my own codes
in many places.
Usage
vcat(vec, sep =" ")
ccat(...)
Arguments
vec, ...
any vector and other arguments to be pasted to together.
sep
the separator to use, see the Details section.
Details
The functions are really just defined as
vcat := function(vec, sep = " ") paste(vec, collapse = sep)
The main motivation for this function has been the easy construction
of a “full GAM formula” from something as simple as
Y ~ ..
The potential use is slightly more general.
Usage
wrapFormula(f, data, wrapString ="s(*)")
Arguments
f
the initial formula; typically something like
Y ~ ..
character string, containing
"*", specifying the wrapping expression to use.
Value
a formula very similar to f; just replacing each
additive term by its wrapped version.
Note
There are limits for this to work correctly; notably the right hand
side of the formula f should not be nested or otherwise
complicated, rather typically just . as in the examples.
Given smoother data (xi,yi) and maybe weights wi,
with multiple xi, use the unique x values, replacing the
y's by their (weighted) mean and updating the weights
accordingly.
Usage
xy.unique.x(x, y, w, fun.mean = mean,...)
Arguments
x, y
numeric vectors of same length. Alternatively, x
can be a ‘xy’ like structure, see xy.coords.
w
numeric vector of non-negative weights – or missing which
corresponds to all weights equal.
,
Werner Stahel [ctb] (Functions: compresid2way(), f.robftest(), last(),
p.scales(), p.dnorm()),
Andreas Ruckstuhl [ctb] (Functions: p.arrows(), p.profileTraces(),
p.res.2x()),
Christian Keller [ctb] (Functions: histBxp(), p.tachoPlot()),
Kjetil Halvorsen [ctb] (Functions: KSd(), ecdf.ksCI()),
Alain Hauser [ctb] (Functions: cairoSwd(), is.whole(),
toLatex.numeric()*),
Christoph Buser [ctb] (to function Duplicated()),
Lorenz Gygax [ctb] (to function p.res.2fact()),
Bill Venables [ctb] (Functions: empty.dimnames(), primes()),
Tony Plate [ctb] (to inv.seq()),
Isabelle Flückiger [ctb],
Marcel Wolbers [ctb],
Markus Keller [ctb],
Sandrine Dudoit [ctb],
Jane Fridlyand [ctb],
Greg Snow [ctb] (to loessDemo()),
Henrik Aa. Nielsen [ctb] (to loessDemo()),
Vincent Carey [ctb],
Ben Bolker [ctb],
Philippe Grosjean [ctb],
Frédéric Ibanez [ctb],
Caterina Savi [ctb],
Charles Geyer [ctb],
Jens Oehlschlägel [ctb]