Статьи

От случайного генератора к случайной функции

В этот уик-энд я написал  пост,  поскольку у меня возникли проблемы с генерацией случайной выборки с использованием R, чтобы воспроизвести одну, полученную соавтором, с SAS (сгенерированным с использованием  Fishman и Moore (1982),  использованным в функции RANUNI). Мне повезло, так как другой автор для этой книги,  Кристроф Дутанг , получил ответ на последний вопрос, который я задал: возможно ли воспроизвести генератор случайных чисел? Да мы можем. И это довольно просто, если вы используете соответствующую библиотеку и соответствующую функцию.

> library(randtoolbox)
> a <- 397204094
> b <- 0
> m <- 2^(31)-1
> set.generator(name="congruRand", mod=m, mult=a, incr=b, seed=123) 
> runif(10)
 [1] 0.7503961 0.3209120 0.1783896 0.9060334 0.3571171
 [6] 0.2211140 0.7864383 0.3980819 0.1246652 0.1876858

Если вы проверите в предыдущем  посте,  это именно то, что дал нам SAS (и что я не мог воспроизвести сам). Но это была только одна часть моих проблем, поскольку цель состояла в том, чтобы на самом деле воспроизвести индексы для обучающей подвыборки по вопросам кредитного скоринга.

Я должен признать, что я никогда не думал об этом раньше: как мы должны написать пример функции? Если значения можно заменить, это нормально, мы просто должны правильно разделить интервал. подобно

> set.seed(95)
> (U=runif(10))
 [1] 0.15171155 0.57584087 0.05309844 0.07044485 0.48887914 0.15276707
 [7] 0.37405684 0.30006096 0.96997126 0.30373498
> set.seed(95)
> (R=sample(0:99,size=10,replace=TRUE))
 [1] 15 57  5  7 48 15 37 30 96 30

Здесь мы просто усекаем значения, полученные из генератора случайных чисел. И это просто отлично. Но как нам написать код для примера  без замен ? Я имею в виду, как вы получаете это:

> (S=sample(0:99,size=10,replace=FALSE))
 [1] 15 57  5  6 46 14 35 27 89 92

Моей первоначальной идеей было использовать следующую технику. Первое значение легко получить: мы просто разбиваем единичный интервал на 100 подразделений (как и раньше, поскольку для первого значения, замена или нет, нам все равно) и видим, в каком интервале находится случайное значение. И мы удаляем это значение из нашего образца. Затем мы разделяем интервал между единицами на 99 и видим, в каком интервале находится случайное значение. Это 10-й? Хорошо, тогда наше второе значение является десятым от нашего образца (первое значение было удалено). Затем мы разделяем интервал между единицами на 98 и т. Д. Код, который я написал для создания этого алгоритма, был следующим:

> mysample1=function(N,unif){
+  n=length(N)
+  size=length(unif)
+  V0=N[trunc(unif[1]*n)+1]
+  N=N[-which(N==V0)]
+  V=V0
+  for(i in 2:length(unif)){
+    V0=N[trunc(unif[i]*(n-i+1))+1]
+    N=N[-which(N==V0)]
+    V=c(V,V0)}
+ return(V)}

К сожалению, я не смог воспроизвести образец, полученный с помощью функции R,

> mysample1(0:99,unif=U)
 [1] 15 58  5  7 49 17 39 31 97 32
> S
 [1] 15 57  5  6 46 14 35 27 89 92

Поскольку  Кристроф  является экспертом по случайным генераторам, я спросил его еще раз. И он придумал следующий код,

> mysample2=function(N,unif){
+   integerset=1:length(N)
+   result=rep(NA,length(unif))
+   for(i in 1:length(unif)){
+     intchosen=integerset[ceiling(U[i]*(length(N)-i+1))]
+     integerset[intchosen]=length(integerset)
+     integerset=integerset[-length(integerset)]
+     result[i]=intchosen}
+   return(N[result])}

который работает просто отлично!

> mysample2(0:99,unif=U)
 [1] 15 57  5  6 46 14 35 27 89 92
> S
 [1] 15 57  5  6 46 14 35 27 89 92

Так что теперь мы можем не только воспроизводить случайные числа, полученные с помощью другого программного обеспечения, мы также можем получать те же индексы выборок с заменой или без нее! Спасибо, Кристоф!

[15 мая]  Обратите внимание, что на самом деле это генератор, используемый в SAS. Чтобы воспроизвести функцию выборки SAS, алгоритм является намного более простым, явно не таким эффективным, поскольку мы генерируем случайную выборку размером 100 (если у нас  есть  100 наблюдений), а затем сохраняем значения, связанные с индексами. из 10 самых маленьких (если мы хотим образец размером 10). Код может быть

> mysample3=function(N,unif,size){
+ q=sort(unif)[size]
+ return(N[U<=q])}

> library(randtoolbox)
> a <- 397204094
> b <- 0
> m <- 2^(31)-1
> set.generator(name="congruRand", mod=m, mult=a, incr=b, seed=123) #OK

> U=runif(100)
> mysample3(1:100,U,size=10)
 [1] 27 37 47 59 60 71 75 82 84 87

Спасибо  Жан-Филиппу  за идею (которая работает).