--Multipart_Mon_Nov_24_14:51:09_1997-1
Content-Type: text/plain; charset=US-ASCII
>>>>> "PaulG" == Paul Gilbert
<pgilbert@bank-banque-canada.ca> writes:
MM> The code is basically in V&R 1 and 2; V&R2 on p.167. I have
it as a
MM> C function that I used to dyn.load into S-plus in order
MM> to prove that S-plus was using it.
PaulG> Could you send this to me. I'd like to try to reproduce some S
PaulG> results in R.
Okay here is a little shar file which contains a little more;
it also includes (the S code for) my experiments on what S-plus's
set.seed(.)
does. The C code must be compiled (cc and ld) ``a la R'' before
dyn.load.
I did not yet test this in R,
but I CC: it to R-devel since others may want to comment.
[appended at the end !]
--------------------------------------------------------
PaulG> Maybe there is some subtlty here that I don't understand, but is
it not
PaulG> possible to just add an optional argument to the functions which
call
PaulG> the random number generator, with the default being the current
PaulG> generator?
Reasons why I wouldn't like this:
- There are too many of these functions; some are quite nested;
- adding an extra argument to them is ugly and not backwards (R-minus)
compatible;
- you shouldn't have to change the R code of your simulation, but rather
change a global option in order to use a different RNG.
- You shouldn't accidentally mix two different generators, since the
resulting RNG has unknown properties.be allowed t
- It isn't too hard, to implement the more flexible scheme,
in the current R implementation.
>From the UI, I'd propose
to use a function instead of .Random.seed,
set.random(generator = "Wichmann-Hill", seed = NULL)
The function (also when called with no arguments) returns a list with the
two components
$generator : character(1) giving the name of the RNG.
$seed : integer vector
AND also accepts such a list as input (alternatively), i.e.
I could
first)
.R.seed <- set.random()
## ....do simulations
later)
set.random(.R.seed)
## ....Re-do simulations indentically
-------------------------------------------------------------------------------
--Multipart_Mon_Nov_24_14:51:09_1997-1
Content-Type: application/octet-stream
Content-Disposition: attachment; filename="runif.shar"
Content-Transfer-Encoding: 7bit
#!/bin/sh
# This is a shell archive (produced by GNU sharutils 4.2).
# To extract the files from this archive, save it to some FILE, remove
# everything before the `!/bin/sh' line above, then type `sh FILE'.
#
# Made on 1997-11-24 14:21 MET by <maechler@sophie>.
# Source directory was `/users/home2/staff/maechler/S/MISC'.
#
# Existing files will *not* be overwritten unless `-c' is specified.
#
# This shar contains:
# length mode name
# ------ ---------- ------------------------------------------
# 7088 -rw-r--r-- runif-tst.S
# 3646 -rw-r--r-- setseed-ex.S
# 1132 -rw-r--r-- ../C-Progs/programs/Runif.c
#
save_IFS="${IFS}"
IFS="${IFS}:"
gettext_dir=FAILED
locale_dir=FAILED
first_param="$1"
for dir in $PATH
do
if test "$gettext_dir" = FAILED && test -f $dir/gettext \
&& ($dir/gettext --version >/dev/null 2>&1)
then
set `$dir/gettext --version 2>&1`
if test "$3" = GNU
then
gettext_dir=$dir
fi
fi
if test "$locale_dir" = FAILED && test -f $dir/shar \
&& ($dir/shar --print-text-domain-dir >/dev/null 2>&1)
then
locale_dir=`$dir/shar --print-text-domain-dir`
fi
done
IFS="$save_IFS"
if test "$locale_dir" = FAILED || test "$gettext_dir" =
FAILED
then
echo=echo
else
TEXTDOMAINDIR=$locale_dir
export TEXTDOMAINDIR
TEXTDOMAIN=sharutils
export TEXTDOMAIN
echo="$gettext_dir/gettext -s"
fi
touch -am 1231235999 $$.touch >/dev/null 2>&1
if test ! -f 1231235999 && test -f $$.touch; then
shar_touch=touch
else
shar_touch=:
echo
$echo 'WARNING: not restoring timestamps. Consider getting and'
$echo "installing GNU \`touch', distributed in GNU File
Utilities..."
echo
fi
rm -f 1231235999 $$.touch
#
if mkdir _sh01181; then
$echo 'x -' 'creating lock directory'
else
$echo 'failed to create lock directory'
exit 1
fi
# ============= runif-tst.S =============if test -f 'runif-tst.S'
&& test "$first_param" != -c; then
$echo 'x -' SKIPPING 'runif-tst.S' '(file already
exists)'
else
$echo 'x -' extracting 'runif-tst.S' '(text)'
sed 's/^X//' << 'SHAR_EOF' > 'runif-tst.S'
&&
##SunOS4: print(dyn.load2("/u/maechler/C/C-for-S/Runif.o"))
print(dyn.load("/u/maechler/C/C-for-S/Runif.o"))
X.C("Init_seed",
X Cval = as.integer(12345),
X Tval = as.integer(987654321))
str(.C("Get_seed", Cval = integer(1), Tval = integer(1)))
X
rU_ numeric(10)
for(i in 1:10) rU[i] <- unlist(.C("Runif", numeric(1)))
rU
str(.C("Get_seed", Cval = integer(1), Tval = integer(1)))
X
Set.seed <- function(Cval = 12345, Tval =987654321)
{
X ## Purpose: Initialize my implementation of the 'S' RNG (following
Ripley)
X ## -------------------------------------------------------------------------
X ## Arguments: Congrval and Tausval, the 2 long integers [== seed]
X ## -------------------------------------------------------------------------
X ## Author: Martin Maechler, Date: 9 May 94, 12:04
X if(!is.loaded(C.symbol("Init_seed")))
X dyn.load("/u/maechler/C/C-for-S/Runif.o")
X invisible(.C("Init_seed",
X Cval = as.integer(Cval),
X Tval = as.integer(Tval)))
}
X
Runif <- function()
{
X ## Purpose: test my runif(.)
X ## -------------------------------------------------------------------------
X ## Author: Martin Maechler, Date: 9 May 94, 11:54
X if(!is.loaded(C.symbol("Runif")))
X dyn.load("/u/maechler/C/C-for-S/Runif.o")
X unlist(.C("Runif", numeric(1)))
}
X
X
tst.runif <- function(Cval = 12345, Tval =987654321)
{
X ## Purpose: test my (following Ripley) implementation of S RNG
X ## return 'next' seeds
X ## -------------------------------------------------------------------------
X ## Arguments: Congrval and Tausval [== seed]
X ## -------------------------------------------------------------------------
X ## Author: Martin Maechler, Date: 9 May 94, 12:04
X if(!is.loaded(C.symbol("Runif")))
X dyn.load("/u/maechler/C/C-for-S/Runif.o")
X .C("Init_seed",
X Cval = as.integer(Cval),
X Tval = as.integer(Tval))
X .C("Runif", numeric(1))
X unlist(.C("Get_seed", Cval = integer(1), Tval = integer(1)))
}
X
###
### Test the 'standard' .Random.seed --- compare with my version
###
X
tst.R.S <- function(nn)
{
X ## Purpose: return many .Random.seed s
X ## -------------------------------------------------------------------------
X ## Arguments: number
X ## -------------------------------------------------------------------------
X ## Author: Martin Maechler, Date: 9 May 94, 11:07
X
X RR <- matrix(0,nrow = 12, ncol= nn)
X for(i in 1:nn) { RR[,i] _ .Random.seed; runif(1) }
X t(RR)
}
X
RR <- tst.R.S(2000)
summary(as.data.frame(RR))#-- Random.seed components
RRt <- apply(RR, 2, table)
RRt2 <- lapply(RRt, function(x) as.integer(names(x)))
str(RRt2, vec.len = 20)
for(j in 1:12) {
X cat("\n j=",j, " 0:63 = R..[[j]]: ", all.equal(0:63,
RRt2[[j]]))
}
X
X
X
X
X
X
##--> only 1, 6, 12 are NOT == 0...63
## [[1]] = 1 5 9 .. 61 (( = 1 + 4*j , j= 0..15 ))
## [[6]] = [[12]] = 0..3
X
X
RS <- rep(0,12)
#.Random.seed _ RS
##>> Error: Invalid random number seed: congruential part must be odd
RS[1] _ 1
X.Random.seed _ RS #--> 1st part == congruential part !
for(i in 1:20){ print(.Random.seed); runif(1)} #- the 2nd part STAYS 0 !!
RS2 <- RS
RS2[12] _ 1
X.Random.seed _ RS2
for(i in 1:20){ print(.Random.seed); runif(1)} #- 1st part == 1st part above!
X
###==> first 6 numbers are one long integer, == congrval (see above)
### second 6 numbers are 2nd long integer == tausval
X
### congrval is made up from 4 + 4* 6 + 3 = 31 bits
### - actually CAN use ALL odd 1:63 for [1] -> 32
bits
### tausval is made up from 5* 6 + 3 = 33 bits (?)
X
X
### The following PROOFS that Runif(.) returns the SAME as runif(.)
### AND that .Random.seed =^= ( Congr.val, Taus.val )
X
if(!exists("polyn.eval"))
X polyn.eval <- function(coef, x)
X {
X ## Purpose: Compute coef[1] + coef[2]*x + ... + coef[p+1]* x^p
X ## x can be any array; result if of same dimensions !
X ## Author: Martin Maechler
X lc <- length(coef)
X if(lc == 0) 0
X else {
X r <- coef[lc]
X if(lc > 1) for(i in (lc - 1):1) r <- coef[i] + r * x
X r
X }
X }
X
qR <- c( 1, rep(0,5))
X.Random.seed <- Rs <- c(qR, qR) #--- 1 0 0 0 0 0 1 0 0 0 0 0
Set.seed(1,1)
for (i in 1:100) {
X CT <- unlist(.C("Get_seed", Cval = integer(1), Tval =
integer(1)))
X RS.C <- .Random.seed[1:6]
X RS.T <- .Random.seed[1:6 +6]
X cat(i,": ")
X ##print(c(CT, RS.C = RS.C, RS.T = RS.T))
X if(CT["Cval"]<0) { CT["Cval"] <- 2^32 +
CT["Cval"]
X cat("'Cval' corrected by adding 2^32\n") }
X if(CT["Tval"]<0) { CT["Tval"] <- 2^32 +
CT["Tval"]
X cat("'Tval' corrected by adding 2^32\n") }
X if(CT["Cval"] != polyn.eval(RS.C, 2^6))
stop("'Congr.val' do not match!")
X if(CT["Tval"] != polyn.eval(RS.T, 2^6)) stop("
'Taus.val' do not match!")
X u1 <- runif(1)
X u2 <- Runif()
X if((u1-u2)!=0) stop("u1 != u2")
}
X
X
X.Random.seed <- Rs <- c(qR, qR) #--- 1 0 0 0 0 0 1 0 0 0 0 0
Set.seed(1,1)
X
runif(1)
##> [1] 4.659873e-05
Runif()
##> [1] 4.659873e-05
unlist(.C("Get_seed", Cval = integer(1), Tval = integer(1)))
##> Cval Tval
##> 69069 131073
X.Random.seed
##> [1] 13 55 16 0 0 0 1 0 32 0 0 0
X
Runif()
##> [1] 0.1106027
runif(1)
##> [1] 0.1106027
unlist(.C("Get_seed", Cval = integer(1), Tval = integer(1)))
##> Cval Tval
##> 475559465 524293
X.Random.seed
##> [1] 41 24 7 22 28 0 5 0 0 2 0 0
X
X
### --- finding out things about set.seed(.) -----
X
for(i in 0:1100){
X cat("set.seed(i=",formatC(i,dig=4),") -> .Random.s =
"); set.seed(i)
X cat(format(.Random.seed),"\n")
}
##- set.seed(i= 0 ) -> .Random.s = 21 14 49 32 43 1 32 22 36 23 28
3
##- set.seed(i= 1 ) -> .Random.s = 21 14 49 48 24 1 32 22 36 23 28
3
##- set.seed(i= 2 ) -> .Random.s = 21 14 49 0 6 1 32 22 36 23 28
3
##- set.seed(i= 3 ) -> .Random.s = 21 14 49 16 51 0 32 22 36 23 28
3
##- set.seed(i= 4 ) -> .Random.s = 21 14 49 32 32 0 32 22 36 23 28
3
##- set.seed(i= 5 ) -> .Random.s = 21 14 49 48 13 0 32 22 36 23 28
3
##- set.seed(i= 6 ) -> .Random.s = 21 14 49 0 59 3 32 22 36 23 28
3
##- set.seed(i= 7 ) -> .Random.s = 21 14 49 16 40 3 32 22 36 23 28
3
##- .........
##- [ 1 2 3 4 5 6 7 8 9 10 12 12 ]
###--> only [4:6] really change !!
###- They have 6 + 6 + 2 = 14 bits each (the ranges are 0:63, 0:63, and 0:3),
###- i.e., 2^14 = 16384 different possibilities.
###- ~~~~ ~~~~~
###===> Period of set.seed(.) =: P | 2^14
X
## test this period
tst.sseed <- function(i)
{
X ## Purpose: test set.seed(.)
X ## -------------------------------------------------------------------------
X ## Author: Martin Maechler, Date: 9 May 94, 11:54
X set.seed(i)
X .Random.seed
}
X
for (i in c(0:7, 11, 17, 20, 23, 40, 53, 59, 67, 76, 97)) {
X cat("i=",i,": ")
X print(all.equal(tst.sseed(i), tst.sseed(i + 2^ 14)))
}
##--> all 'TRUE'
##
## same thing for 2^13, 2^12, 2^11, 2^10 ! --
##
## Test a full period for 2^10 = 1024 :
X
for (i in 0:(2^10 - 1))
X if(!is.all.equal(tst.sseed(i), tst.sseed(i + 2^ 10)))
X cat("i=",i," DIFFERENT ! ")
## -> nothing (i.e., all equal !)
X
###==> Full period is from set.seed(0) to set.seed(1023) ! --
SHAR_EOF
$shar_touch -am 1124142097 'runif-tst.S' &&
chmod 0644 'runif-tst.S' ||
$echo 'restore of' 'runif-tst.S' 'failed'
if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' )
>/dev/null 2>&1 \
&& ( md5sum --version 2>&1 | grep -v 'textutils 1.12' )
>/dev/null; then
md5sum -c << SHAR_EOF >/dev/null 2>&1 \
|| $echo 'runif-tst.S:' 'MD5 check failed'
35e91fce41f56de76bab514ecfc5c840 runif-tst.S
SHAR_EOF
else
shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c <
'runif-tst.S'`"
test 7088 -eq "$shar_count" ||
$echo 'runif-tst.S:' 'original size' '7088,'
'current size' "$shar_count!"
fi
fi
# ============= setseed-ex.S =============if test -f 'setseed-ex.S'
&& test "$first_param" != -c; then
$echo 'x -' SKIPPING 'setseed-ex.S' '(file already
exists)'
else
$echo 'x -' extracting 'setseed-ex.S' '(text)'
sed 's/^X//' << 'SHAR_EOF' > 'setseed-ex.S'
&&
###--> See also (the much OLDER) ./runif-tst.S
### ~~~~~~~~~~~~~
X
ss <- function(i) { set.seed(i); .Random.seed }
s0 <- ss(0)
for(i in -1000:4000)
X if(any(ss(i)[-(4:6)] != s0[-(4:6)])) cat("i=",i,"
.R.s=",.Random.seed,"\n")
##--> no output, i.e., it's really
## only .Random.seed[c(4,5,6)] which are modified by set.seed(.)
X
ni <- length(i.set <- -2:2000)
m456 <- matrix(NA, ni, 3, dimnames=list(paste(i.set),paste(4:6)))
for(i in i.set) m456[paste(i),] <- ss(i)[4:6]
t(m456[1:25,])
X
apply(m456, 2, function(x) sort(unique(x)))
##- $"4":
##- [1] 0 16 32 48
##-
##- $"5":
##-[1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23
##[25] 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
##[49] 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
##-
##- $"6":
##- [1] 0 1 2 3
len456 <- apply(m456, 2, function(x) length(unique(x)))
len456
prod(len456) #-- 1024 : no more than these many different
X
d456 <- dim(table(m456[,"4"], m456[,"5"],
m456[,"6"]))#- and these all appear ..
d456
d456 == len456 # Yes !
X
all(m456[,"4"] == rep(16*(0:3), length=ni))#-> TRUE
X
for(i in -4:10) print(ss(i)[4] == 16*(i+2)%%4)
X
table(m456[,"6"]) #-- only 0, 1, 2, 3, but not as regularly ...
##- 0 1 2 3
##- 502 501 500 500
X
##---- mI contains .Random.seed[4:6] after set.seed(i = 1:1024) : ----
mI <- m456[paste(1:1024),]
dim(table(mI[,"4"], mI[,"5"], mI[,"6"]))#- and
these all appear ..
X
all(ss(0) == ss(1024)) #-- TRUE
for(i in -100:2000)
X if(any(ss(i) != ss(i + 1024))) cat("i=",i,"
.R.s=",.Random.seed,"\n")
#- no output, i.e., it's true
X
# the part of the 'congrval'
all(table(mI %*% 2^(6*0:2)) == 1)# TRUE
X
X
## Return 'congrval' part of .Random.seed after set.seed(.):
## need 'polyn.eval' from my goodies:
ssC <- function(i){ set.seed(i);polyn.eval(.Random.seed[ 1:6],2^6)}
ssT <- function(i){ set.seed(i);polyn.eval(.Random.seed[6+1:6],2^6)}
X
##-- All the 'congrval' that set.seed(.) produces ..
Cs <- ii <- 1:1024
for(i in ii) Cs[i] <- ssC(i)
dim(tCs <- table(Cs))# = 1024
all(tCs == 1) # TRUE (each once)
X
summary(Cs / 2^32)
##- Min. 1st Qu. Median Mean 3rd Qu. Max.
##- 4.694e-05 0.2498 0.4996 0.4996 0.7493 0.9991
X
dCs <- unique(diff(srt.C <- sort(Cs)))
dCs##-> 4194304 only one value !
log(dCs, base =2)# == 22
X
C1 <- srt.C[1]
all(srt.C == C1 + (ii-1) * 2^22) #- TRUE !
X
## the values are not quite symmetrical in {0,..,2^32-1} :
2^32 - srt.C[1024]# 3992683
C1 # 201621
X
###------------ Martin M鋍hler's Summary : ------------------------
X
###--- Mon, Sept. 15, 1997 ----------
X
##- Unfortunately, we don't know what exactly happens in
##- the S-plus code for set.seed(i), namely, .C("setseed",
as.integer(i))
##- However, I have found out the following which should suffice :
##-
##- 1) Only .Random.seed[4:6] are varyied by set.seed, i.e.,
##- after set.seed(i), .Random.seed always is
##-
##- 1 2 3 4 5 6 7 8 9 10 11 12
##- [1] 21 14 49 * * * 32 22 36 23 28 3
##- X4 X5 X6
##-
##- which means that only the 'congrval' {Venables & Ripley, 2nd
ed., p.166-167}
##- is varied
##-
##- 2) X4 is always == 16*((i+2) %% 4)
##- X5 takes (all) values in 0:63 [also with some periodicity]
##- X6 takes the values of 0:3
##-
##- 3) set.seed(i) <==> set.seed(i + n*1024) for all integers n, i
##-
##- 4) { set.seed(i) ; i \in {integers}}
##- results in a set of 1024 different values of 'congrval'.
##- These values are EXACTLY EQUISPACED lying in {0,..,2^32}.
##-
##- If C[1:1024] are these sorted values, we have
##-
##- C[i] == 201621 + (i-1)* 2^22 , for i in 1:1024
##-
X
SHAR_EOF
$shar_touch -am 1124095797 'setseed-ex.S' &&
chmod 0644 'setseed-ex.S' ||
$echo 'restore of' 'setseed-ex.S' 'failed'
if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' )
>/dev/null 2>&1 \
&& ( md5sum --version 2>&1 | grep -v 'textutils 1.12' )
>/dev/null; then
md5sum -c << SHAR_EOF >/dev/null 2>&1 \
|| $echo 'setseed-ex.S:' 'MD5 check failed'
c2d03d4318c03a8f2ed8b4f0b656c647 setseed-ex.S
SHAR_EOF
else
shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c <
'setseed-ex.S'`"
test 3646 -eq "$shar_count" ||
$echo 'setseed-ex.S:' 'original size' '3646,'
'current size' "$shar_count!"
fi
fi
# ============= ../C-Progs/programs/Runif.c =============if test ! -d
'..'; then
$echo 'x -' 'creating directory' '..'
mkdir '..'
fi
if test ! -d '../C-Progs'; then
$echo 'x -' 'creating directory' '../C-Progs'
mkdir '../C-Progs'
fi
if test ! -d '../C-Progs/programs'; then
$echo 'x -' 'creating directory' '../C-Progs/programs'
mkdir '../C-Progs/programs'
fi
if test -f '../C-Progs/programs/Runif.c' && test
"$first_param" != -c; then
$echo 'x -' SKIPPING '../C-Progs/programs/Runif.c' '(file
already exists)'
else
$echo 'x -' extracting '../C-Progs/programs/Runif.c'
'(text)'
sed 's/^X//' << 'SHAR_EOF' >
'../C-Progs/programs/Runif.c' &&
/*--- this code is from a posting of
X Brian Ripley to the S-news, 7 Aug 92,
X and a slight correction from Allen Mcintosh
<mcintosh@larch.bellcore.com>
X + a correction of mine (B.R. wrote '69096' instead of
'69069'))
X
X --> /u/maechler/tex/o.peopl/Ripley-S-RNG.tex
*/
X
/*--- This code is perfectly tested out via S-code
X
X in --> /u/maechler/S/MISC/runif-tst.S
*/
X
unsigned long congr_MM, taus_MM;
X
static double xuni()
{
X unsigned long n, lambda = 69069;
X congr_MM *= lambda;
X taus_MM ^= taus_MM >> 15;
X taus_MM ^= taus_MM << 17;
/* Orig (Splus <= 3.2):
X n = taus_MM ^ congr_MM;
X return ( (n>>1) / 2147483648.);
X * New: */
X n = (taus_MM ^ congr_MM) >> 1;
X if(n) { /* n != 0 */ return( n / 2147483648.);
X } else return(xuni());
}
X
X
/*--- The following are to be dyn.load(..)ed and called from S(plus) --*/
X
void Runif (double* res) { *res = xuni(); }
X
void Init_seed (long* Cval, long* Tval)
{
X congr_MM = (unsigned long)*Cval;
X taus_MM = (unsigned long)*Tval;
}
X
void Get_seed (long* Cval, long* Tval)
{
X *Cval = (long) congr_MM;
X *Tval = (long) taus_MM;
}
SHAR_EOF
$shar_touch -am 1124141497 '../C-Progs/programs/Runif.c' &&
chmod 0644 '../C-Progs/programs/Runif.c' ||
$echo 'restore of' '../C-Progs/programs/Runif.c'
'failed'
if ( md5sum --help 2>&1 | grep 'sage: md5sum \[' )
>/dev/null 2>&1 \
&& ( md5sum --version 2>&1 | grep -v 'textutils 1.12' )
>/dev/null; then
md5sum -c << SHAR_EOF >/dev/null 2>&1 \
|| $echo '../C-Progs/programs/Runif.c:' 'MD5 check failed'
f5888c70e55b647a34d8bee3aae3b4ad ../C-Progs/programs/Runif.c
SHAR_EOF
else
shar_count="`LC_ALL= LC_CTYPE= LANG= wc -c <
'../C-Progs/programs/Runif.c'`"
test 1132 -eq "$shar_count" ||
$echo '../C-Progs/programs/Runif.c:' 'original size'
'1132,' 'current size' "$shar_count!"
fi
fi
rm -fr _sh01181
exit 0
--Multipart_Mon_Nov_24_14:51:09_1997-1
Content-Type: text/plain; charset=US-ASCII
--Multipart_Mon_Nov_24_14:51:09_1997-1--
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=---
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To:
r-devel-request@stat.math.ethz.ch
=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=---
