RNG in R and Splus 3.4
Prof. Ripley asked the details of the example.
We were doing parametric bootstrap, so it is similar to simulation.
Anyway here is the details.
We start with a sample of 19 positive numbers. We know the sample
is from truncated exp(0.3)...only the truncation point, theta, is unknown.
In other words, the sample can be generated from something like
x1 <- rexp(100, rate=0.3)
x2 <- x1[x1<theta] # the theta are big enough so that length(x2)>19
x3 <- x2[1:19]
And the true theta was 4, hidden from the students.
The estimate of theta from that particular sample was 3.87058501,
the max of the sample.
Next we wanted to estimate the bias of the max, using bootstrap.
We did this:
y1<-rexp(4000, rate=0.3)
y2<- y1[y1<3.87058501]
length(y2) # to make sure y2 is longer than 1900 or whatever..so we are not
# recycling
y3<-y2[1:1900]
# now y3 contains 100 bootstrap samples of 19 each
y4 <- matrix(y3, ncol=19, byrow=T) # actually some of students did
# y4 <- t( matrix(y3, ncol=100) )
maxboot <- apply(y4, 1, max)
biasboot <- mean(maxboot) - 3.87058501
We did this for more than 100 runs, depend how much memory you have on your
machine, you can try put a few zeros to the above code. but you get the idea.
The final result, biasboot, is very different when using R or Splus 3.4
Roughly 2 times (biasboot from R) = (biasboot from Splus 3.4)
In this particula case the bias was 3.87058501-4. R seem to be better
than Splus which give a biasboot that is too small (or too big in abs value)
But more analysis using simulation reveals otherwise.
We then tried the simulation and that was reported earlier....
We use Rjune on the Win95 and Win NT machine.
Mai Z
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
On Mon, 9 Nov 1998, Mai Zhou wrote:> RNG in R and Splus 3.4 > > Prof. Ripley asked the details of the example. > We were doing parametric bootstrap, so it is similar to simulation. > Anyway here is the details.The following code replicates this without needing any specific hand-coding of numbers each time it is run (except for the slight chance that y2 will be shorter than 1900 elements) x1 <- rexp(100, rate=0.3) theta<-4 x2 <- x1[x1<theta] x3 <- x2[1:19] y1<-rexp(4000, rate=0.3) y2<-y1[y1< max(x3)] y3<-y2[1:1900] y4 <- matrix(y3, ncol=19, byrow=T) maxboot <- apply(y4, 1, max) biasboot <- mean(maxboot) - max(x3) biasboot I ran this 1000 times on S-PLUS 3.4 for Solaris, Rpre0.63 for Solaris, Guido's R 0.61.3 for Windows (NT) and R0.61.3 for Linux (the oldest R I have lying around) and got S-PLUS -0.3040023 Rpre0.63 -0.3083820 R 0.61.3 -0.3054932 WinNT -0.307758 I don't have easy access to RJune, but it seems that on the other current and recent Rs the example works consistently. Anyone else tried this on RJune? -thomas Thomas Lumley Assistant Professor, Biostatistics University of Washington, Seattle. -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.- 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 _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Thomas, Interesting, can you try the following too?
b <- function(x, nboot, theta) {
y1<-rexp(2*length(x)*nboot, rate=0.3)
y2<-y1[y1<theta]
y3<-y2[1:(length(x)*nboot)]
y4<- matrix(y3, ncol=length(x), byrow=T)
mean(apply(y4, 1, max))- theta
}
I got arround -0.155 on Rjune.
Where x is a length 19 vector, theta=4
One typical run is > b(rep(1,19), 1500, 4)
[1] -0.1585402
While on Splus I got> b(rep(1,19), 1500, 4)
[1] -0.3516429
Mai Z>
> I ran this 1000 times on S-PLUS 3.4 for Solaris, Rpre0.63 for Solaris,
> Guido's R 0.61.3 for Windows (NT) and R0.61.3 for Linux (the oldest R
I
> have lying around) and got
>
> S-PLUS -0.3040023
> Rpre0.63 -0.3083820
> R 0.61.3 -0.3054932
> WinNT -0.307758
>
> I don't have easy access to RJune, but it seems that on the other
current
> and recent Rs the example works consistently.
>
> Anyone else tried this on RJune?
>
> -thomas
>
>
> Thomas Lumley
> Assistant Professor, Biostatistics
> University of Washington, Seattle.
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help 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-help-request at
stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._