Replicating a random normal generated in SAS (rancor) in R based on the same seed?
To do this, you need two things:
- Seed used to generate a random number
- Formula used to generate a random number
SAS uses for rannor
(and I think also for rand
, but I have not seen confirmation of this), the following algorithm (found in Psuedo-Random Numbers: Out of Uniform , Robert Johnson and Hui Liu):
rannor_result = (−2* log(U1))**.5*cos(2*constant('pi')*U2)
where U1 and U2 are two numbers from a uniform numerical stream. (A second number could also be obtained, but that number is discarded as far as I can tell.)
See the following SAS file:
data test1;
U1 = ranuni(7);
U2 = ranuni(7);
X1 = (−2* log(U1))**.5*cos(2*constant('pi')*U2);
U1 = ranuni(7);
U2 = ranuni(7);
X2 = (−2* log(U1))**.5*cos(2*constant('pi')*U2);
run;
data test2;
x1 = rannor(7);
x2 = rannor(7);
run;
test1
and test2
have the same meaning for the runner numbers.
I doubt R and SAS share common PRNG algorithms, especially if using one is rannor
n't very good (you should use rand
, which is much better if using the Mersenne Twister). However, you can easily query SAS to output the seeds it used - as long as you are outputting RANUNI results.
You can ask SAS to do it like this:
data rands_uni;
seed=7;
do _i_=1 to 10;
seed1 = seed;
call ranuni(seed,x1);
seed2=seed;
call ranuni(seed,x2);
output;
end;
run;
From this, you could calculate the runner result in R (i.e. from x1 and x2). I'm including the seeds there for reference - maybe R has the ability to use them. The article above mentions the SAS algorithm used for ranuni
, but it also claims that you cannot fully reproduce it due to some fixes (or maybe due to floating point issues?).
source to share
I am testing both, it doesn't seem to be the case.
data _null_;
CALL STREAMINIT(1245);
do i=1 to 10;
x2 = rand('NORMAL',0,1);
put x2;
end;
run;
/*
-1.295643934
0.0651085355
-0.391011151
1.1024205822
2.0712099724
-0.296163515
0.898941256
-0.584950901
0.2027164546
-0.986522152*/
set.seed(1245)
df<-rnorm(10,0,1)
df
[1] 0.41914007 -1.45583361 -1.45588008 -1.01564637 -0.53309914 -0.68360608 0.32590880
[8] 0.92256137 -0.05085288 0.29722361
source to share
Proc IML;
Call RandSeed(6);
U1 = J(6,1);
U2 = J(6,1);
/* Weibull Distribution */
Call RandGen(U1,"WEIBULL",1.7737152,13.164695);
Print "Weibull Distribution:" U1;
/* Unifrom Distribution */
Call RandGen(U2,"UNIFORM",0,1);
Print "Uniform Distribution:" U2;
Submit/R;
set.seed(6);
## Weibull Distribution;
U1 <- rweibull(6,1.7737152,13.164695);
cat("Weibull Distribution:", U1, "\n");
## Uniform Distribution;
U2 <- runif(6,0,1);
cat("Unifrom Distribution:", U2, "\n");
EndSubmit;
Quit;
source to share