Re: [Rd] A bug in the R Mersenne Twister (RNG) code?
> Dirk Eddelbuettel > on Wed, 31 Aug 2016 10:30:07 -0500 writes: > On 30 August 2016 at 18:29, Duncan Murdoch wrote: > I don't see evidence of a bug. There have been several > versions of the MT; we may be using a different version > than you are. Ours is the 1999/10/28 version; the web > page you cite uses one from 2002. > > Perhaps the newer version fixes some problems, and then > it would be worth considering a change. But changing > the default RNG definitely introduces problems in > reproducibility, so it's not obvious that we would do it. > Yep. FWIW the GNU GSL adopted the 2002 version a while ago > too. Quoting from > https://www.gnu.org/software/gsl/manual/html_node/Random-number-generator-algorithms.html > Generator: gsl_rng_mt19937 >The MT19937 generator of Makoto Matsumoto and Takuji > Nishimura is a variant of the twisted generalized feedback > shift-register algorithm, and is known as the “Mersenne > Twister” generator. It has a Mersenne prime period of > 2^19937 - 1 (about 10^6000) and is equi-distributed in 623 > dimensions. It has passed the DIEHARD statistical > tests. It uses 624 words of state per generator and is > comparable in speed to the other generators. The original > generator used a default seed of 4357 and choosing s equal > to zero in gsl_rng_set reproduces this. Later versions > switched to 5489 as the default seed, you can choose this > explicitly via gsl_rng_set instead if you require it. >For more information see, > Makoto Matsumoto and Takuji Nishimura, “Mersenne > Twister: A 623-dimensionally equidistributed uniform > pseudorandom number generator”. ACM Transactions on > Modeling and Computer Simulation, Vol. 8, No. 1 > (Jan. 1998), Pages 3–30 The generator gsl_rng_mt19937 uses > the second revision of the seeding procedure published by > the two authors above in 2002. The original seeding > procedures could cause spurious artifacts for some seed > values. They are still available through the alternative > generators gsl_rng_mt19937_1999 and gsl_rng_mt19937_1998. > Note the last sentence here. > This is all somewhat technical code, so when I noticed the > above I could never figure what exactly R was doing in its > implementation. But "innocent until proven guilty" -- a > sufficient number of people ought to have looked at this > -- so I saw no need to pursue this further. > Dirk This interesting thread went to sleep a bit early. Let me summarize my view (another R-core after Duncan's) on this: a) R's reference documentation, easily accessed with one of ?Random, ?RNG, ?.Random.seed (and more) clearly indicates the reference for our "Mersenne-Twister" (MT) as the ' Matsumoto, M. and Nishimura, T. (1998) ... ' publication. and we are providing that. As Duncan said, there is no bug. b) The extra information by Mark and notably Dirk indicates that "the litterature" showed that the 1998 version of MT has rare problems and the GSL, notably uses the 2002 version of MT. c) I think it would be nice if we could provide that as well as another RNGkind(). I for one would be willing to integrate a well written patch proposal (what others call PR for "pull request", and I'd also call "code donation") into the R sources. Well-written for me would include to re-use / share code with the 1998 version as much as possible. d) If we had both kinds, say "Mersenne-Twister" and "Mersenne-Twister_2002", we (the R community, not R core necessarily) could conduct experiments about the differences... and can contempltate the pros and cons of a potential switch of default. e) A bit tangential to this: Nowadays, there also exist faster gamma variate RNGs .. giving different random values - for rgamma() but also several other derived RNGs, notably rchisq(), and even more in other CRAN packages. With code donation there, we could introduce a new argument 'gamma.kind' to the set.seed(see, kind = NULL, normal.kind = NULL) RNGkind ( kind = NULL, normal.kind = NULL) functions. (There, currently I'd guess we would not change the default). --- Note that the above includes my guess that R-core would not take action unless we get patch proposals. (But then, I'm happy if my guess is wrong here ...) Martin Maechler, ETH Zurich __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] The use of match.fun
Dear gurus, I was utterly surprised to learn that one of my examples illustrating the need of match.fun() doesn't give me the expected result. center <- function(x,FUN) FUN(x) center(1:10, mean) mean <- 4 center(1:10, mean) Used to give me the error message "could not find function FUN". Now it just works, even though I didn't expect it to. I believe this is at least partially linked to a change in how R finds functions. Now I'm not sure any more whether match.fun() actually has any use any longer, and if so, in which cases it prevents things going wrong. I've tried to find an example where this went wrong, but couldn't find one. Any pointer to what happened here is greatly appreciated. I've checked the NEWS, but I'm not smart enough to find the relevant bits and piece it together. Thank you in advance Cheers Joris -- Joris Meys Statistical consultant Ghent University Faculty of Bioscience Engineering Department of Mathematical Modelling, Statistics and Bio-Informatics tel : +32 (0)9 264 61 79 joris.m...@ugent.be --- Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php [[alternative HTML version deleted]] __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] The use of match.fun
I am not able to replicate this: > center <- function(x,FUN) FUN(x) > center(1:10, mean) [1] 5.5 > mean <- 4 > center(1:10, mean) Error in center(1:10, mean) : could not find function "FUN" Using a fresh install of version 3.3.1 under MacOS, and tested before with 3.3.0 with the same result. On Tue, Sep 6, 2016 at 4:25 PM, Joris Meys wrote: > Dear gurus, > > I was utterly surprised to learn that one of my examples illustrating the > need of match.fun() doesn't give me the expected result. > > center <- function(x,FUN) FUN(x) > center(1:10, mean) > mean <- 4 > center(1:10, mean) > > Used to give me the error message "could not find function FUN". Now it > just works, even though I didn't expect it to. I believe this is at least > partially linked to a change in how R finds functions. > > Now I'm not sure any more whether match.fun() actually has any use any > longer, and if so, in which cases it prevents things going wrong. > > I've tried to find an example where this went wrong, but couldn't find one. > Any pointer to what happened here is greatly appreciated. I've checked the > NEWS, but I'm not smart enough to find the relevant bits and piece it > together. > > Thank you in advance > Cheers > Joris > > -- > Joris Meys > Statistical consultant > > Ghent University > Faculty of Bioscience Engineering > Department of Mathematical Modelling, Statistics and Bio-Informatics > > tel : +32 (0)9 264 61 79 > joris.m...@ugent.be > --- > Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php > > [[alternative HTML version deleted]] > > __ > R-devel@r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel > -- Adrian Dusa University of Bucharest Romanian Social Data Archive Soseaua Panduri nr.90 050663 Bucharest sector 5 Romania [[alternative HTML version deleted]] __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] The use of match.fun
Sorry, I am being a daft idiot again. Turns out that somehow I had an object FUN in my global environment, which explains why I didn't get the result I expected. So please ignore my question and burn my mail on a stake. Cheers Joris On Tue, Sep 6, 2016 at 3:35 PM, Adrian Dușa wrote: > I am not able to replicate this: > > > center <- function(x,FUN) FUN(x) > > center(1:10, mean) > [1] 5.5 > > mean <- 4 > > center(1:10, mean) > Error in center(1:10, mean) : could not find function "FUN" > > Using a fresh install of version 3.3.1 under MacOS, and tested before with > 3.3.0 with the same result. > > > On Tue, Sep 6, 2016 at 4:25 PM, Joris Meys wrote: > >> Dear gurus, >> >> I was utterly surprised to learn that one of my examples illustrating the >> need of match.fun() doesn't give me the expected result. >> >> center <- function(x,FUN) FUN(x) >> center(1:10, mean) >> mean <- 4 >> center(1:10, mean) >> >> Used to give me the error message "could not find function FUN". Now it >> just works, even though I didn't expect it to. I believe this is at least >> partially linked to a change in how R finds functions. >> >> Now I'm not sure any more whether match.fun() actually has any use any >> longer, and if so, in which cases it prevents things going wrong. >> >> I've tried to find an example where this went wrong, but couldn't find >> one. >> Any pointer to what happened here is greatly appreciated. I've checked the >> NEWS, but I'm not smart enough to find the relevant bits and piece it >> together. >> >> Thank you in advance >> Cheers >> Joris >> >> -- >> Joris Meys >> Statistical consultant >> >> Ghent University >> Faculty of Bioscience Engineering >> Department of Mathematical Modelling, Statistics and Bio-Informatics >> >> tel : +32 (0)9 264 61 79 >> joris.m...@ugent.be >> --- >> Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php >> >> [[alternative HTML version deleted]] >> >> __ >> R-devel@r-project.org mailing list >> https://stat.ethz.ch/mailman/listinfo/r-devel >> > > > > -- > Adrian Dusa > University of Bucharest > Romanian Social Data Archive > Soseaua Panduri nr.90 > 050663 Bucharest sector 5 > Romania > -- Joris Meys Statistical consultant Ghent University Faculty of Bioscience Engineering Department of Mathematical Modelling, Statistics and Bio-Informatics tel : +32 (0)9 264 61 79 joris.m...@ugent.be --- Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php [[alternative HTML version deleted]] __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] R (development) changes in arith, logic, relop with (0-extent) arrays
Yesterday, changes to R's development version were committed, relating
to arithmetic, logic ('&' and '|') and
comparison/relational ('<', '==') binary operators
which in NEWS are described as
SIGNIFICANT USER-VISIBLE CHANGES:
[.]
• Arithmetic, logic (‘&’, ‘|’) and comparison (aka
‘relational’, e.g., ‘<’, ‘==’) operations with arrays now
behave consistently, notably for arrays of length zero.
Arithmetic between length-1 arrays and longer non-arrays had
silently dropped the array attributes and recycled. This
now gives a warning and will signal an error in the future,
as it has always for logic and comparison operations in
these cases (e.g., compare ‘matrix(1,1) + 2:3’ and
‘matrix(1,1) < 2:3’).
As the above "visually suggests" one could think of the changes
falling mainly two groups,
1) <0-extent array> (op)
2) <1-extent array> (arith)
These changes are partly non-back compatible and may break
existing code. We believe that the internal consistency gained
from the changes is worth the few places with problems.
We expect some package maintainers (10-20, or even more?) need
to adapt their code.
Case '2)' above mainly results in a new warning, e.g.,
> matrix(1,1) + 1:2
[1] 2 3
Warning message:
In matrix(1, 1) + 1:2 :
dropping dim() of array of length one. Will become ERROR
>
whereas '1)' gives errors in cases the result silently was a
vector of length zero, or also keeps array (dim & dimnames) in
cases these were silently dropped.
The following is a "heavily" commented R script showing (all ?)
the important cases with changes :
(m <- cbind(a=1[0], b=2[0]))
Lm <- m; storage.mode(Lm) <- "logical"
Im <- m; storage.mode(Im) <- "integer"
## 1. -
try( m & NULL ) # in R <= 3.3.x :
## Error in m & NULL :
## operations are possible only for numeric, logical or complex types
##
## gives 'Lm' in R >= 3.4.0
## 2. -
m + 2:3 ## gave numeric(0), now remains matrix identical to m
Im + 2:3 ## gave integer(0), now remains matrix identical to Im (integer)
m > 1 ## gave logical(0), now remains matrix identical to Lm (logical)
m > 0.1[0] ## ditto
m > NULL ## ditto
## 3. -
mm <- m[,c(1:2,2:1,2)]
try( m == mm ) ## now gives error "non-conformable arrays",
## but gave logical(0) in R <= 3.3.x
## 4. -
str( Im + NULL) ## gave "num", now gives "int"
## 5. -
## special case for arithmetic w/ length-1 array
(m1 <- matrix(1,1,1, dimnames=list("Ro","col")))
(m2 <- matrix(1,2,1, dimnames=list(c("A","B"),"col")))
m1 + 1:2 # -> 2:3 but now with warning to "become ERROR"
tools::assertError(m1 & 1:2)# ERR: dims [product 1] do not match the length of
object [2]
tools::assertError(m1 < 1:2)# ERR: (ditto)
##
## non-0-length arrays combined with {NULL or double() or ...} *fail*
### Length-1 arrays: Arithmetic with |vectors| > 1 treated array as scalar
m1 + NULL # gave numeric(0) in R <= 3.3.x --- still, *but* w/ warning to "be
ERROR"
try(m1 > NULL)# gave logical(0) in R <= 3.3.x --- an *error* now in R >=
3.4.0
tools::assertError(m1 & NULL)# gave and gives error
tools::assertError(m1 | double())# ditto
## m2 was slightly different:
tools::assertError(m2 + NULL)
tools::assertError(m2 & NULL)
try(m2 == NULL) ## was logical(0) in R <= 3.3.x; now error as above!
Note that in R's own 'nls' sources, there was one case of
situation '2)' above, i.e. a 1x1-matrix was used as a "scalar".
In such cases, you should explicitly coerce it to a vector,
either ("self-explainingly") by as.vector(.), or as I did in
the nls case by c(.) : The latter is much less
self-explaining, but nicer to read in mathematical formulae, and
currently also more efficient because it is a .Primitive.
Please use R-devel with your code, and let us know if you see
effects that seem adverse.
In some case where R-devel now gives an error but did not
previously, we could contemplate giving another "warning
'to become ERROR'" if there was too much breakage, though
I don't expect that.
For the R Core Team,
Martin Maechler,
ETH Zurich
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
