[R] How to properly re-set a saved seed? I've got the answer, but no explanation

[Paul Johnson]

Hello, happy new year.

I've come back to a problem from last spring. I need to understand
what what is the technically correct method to save a seed and then
re-set it.  Although this example arises in the context of MT19937,
I'm needing to do this with other generators as well, including
L'Ecuyer streams.

The puzzle is this comment in ?Random: "‘set.seed’ is the recommended
way to specify seeds."

What I did not understand before, and can't make sense of now, is that
set.seed does not "re-set" a saved seed.  Here's my working example:


## Paul Johnson
## April 20, 2011

## If you've never dug into the R random number generator and its use of the
## default MT19937, this might be informative.  It will also help you
## avoid a time-consuming mistake that I've made recently.

## I've been developing a simulation and I want to save and restore random
## streams exactly. I got that idea from John Chambers Software for
## Data Analysis and I studied his functions to see how he did it.
## The problem, as we see below, is that set.seed() doesn't do what I expected,
## and I feel lucky to have noticed it now, rather than later.

## I wish set.seed did work the way I want, though, and that's why I'm
## writing here, to see if anybody else wishes the same.

## Here's a puzzle.  Try this:
set.seed(444)
s1 <- .Random.seed
runif(1)
rnorm(1)

set.seed(444)
runif(1)
rnorm(1)
## Those matched. Good

## Re-insert the "saved seed" s1
set.seed(s1)
runif(1)
rnorm(1)

## Why don't the draws match? I "put back" the seed, didn't I?
## Hm. Try again
set.seed(s1)
runif(1)
rnorm(1)

## Why did those match? But neither matches cases 1 and 2.

## I was baffled & discouraged.

## The help page for random numbers says:
##  ‘set.seed’ is the recommended way to specify seeds.

## But, it doesn't say something like

# set.seed(s1)

## is supposed to work. Ah. My mistake.

############ Here's the fix and the puzzle #######
## To re-set the generator to its position at s1, it is required
## to do what the help page seems to say we should not.

.Random.seed <- s1

runif(1)
#########################################

-- 
Paul E. Johnson
Professor, Political Science
1541 Lilac Lane, Room 504
University of Kansas