Воспроизводимость и случайность

Со  Стефаном Туффери мы на этой неделе работали над главой книги под названием «  Статистическое обучение в актуарной науке» . Глава должна быть основана на R-функциях, и мы хотели воспроизвести некоторые результаты, которые он ранее получил с помощью SAS. Хорошая вещь заключается в том, что даже сложные функции (логистическая регрессия, деревья регрессии и т. Д.) Дают одинаковые результаты. Но мы обнаружили проблему, которую не смогли решить: генерирование идентичных обучающих наборов наблюдений … Из 1000 строк мы подобрали около 600 строк. Проблема в том, что мы не могли генерировать одинаковые наборы индексов как с R, так и с SAS. Даже используя аналогичные генераторы случайных чисел … (кроме случаев, когда мы хотим извлечь 1 или 2 строки, не более).

Давайте попробуем объяснить, что происходит (на основе кода, созданного  Стефаном ). Согласно  Eubank (2010) , в SAS есть (как минимум) два генератора случайных чисел,

• Фишман и Мур (1982)  используется для функции RANUNI
• Мерсенн-Твистер,  используемый для функции RAND, основанный на  Мацумото и Нисимуре (1997)

Например, для функции RAND, если мы сгенерируем образец Гаусса — с помощью алгоритма Мерсенна-Твистера — код должен быть

%LET SEED =6;
%LET NREP=10;
DATA TESTRANDOM ;
  CALL STREAMINIT(&SEED);
  DO REP = 1 TO &NREP;
   X = RAND ('NORMAL'); 
   OUTPUT;
  END;
RUN;
PROC PRINT DATA = TESTRANDOM ;
RUN ;

И мы здесь

Если мы хотим единообразный образец, он должен быть

%LET SEED =6;
%LET NREP=10;
DATA TESTRANDM ;
CALL STREAMINIT(&SEED);
DO REP = 1 TO &NREP;
  X = RAND ('UNIFORM'); 
  OUTPUT;
END;
RUN;
PROC PRINT DATA = TESTRANDOM ;
RUN ;

Хорошая новость (пока что) в том, что Mersenne-Twister был закодирован в R, в  функции RNGkind

> RNGkind("Mersenne")
> set.seed(6)
> runif(10)
 [1] 0.64357277 0.91590888 0.09523258 0.29537280
 [5] 0.76993170 0.25589353 0.51789573 0.67784993
 [9] 0.14722782 0.70052604

> a = 397204094      # RANUNI multiplier
> seed = 123         # seed
> n = 10             # sample size
> m = (2^31) - 1     # period
> for (i in (1:n-1)) {
+ seed = (a*seed)%%m
+ u = seed / m
+ print(u)
+ }
[1] 0.7503961
[1] 0.3209121
[1] 0.3453204
[1] 0.2683455
[1] 0.241798
[1] 0.9888181
[1] 0.3943279
[1] 0.9710172
[1] 0.001632214
[1] 0.921537

%LET SEED =123;
%LET NREP=10;
DATA FRANUNI (KEEP = x) ;
seed = &SEED ;
DO REP = 1 TO &NREP;
  CALL RANUNI(seed, x);
  OUTPUT;
END;
RUN;
PROC PRINT DATA = FRANUNI ;
RUN ;