В этот уик-энд я написал пост, поскольку у меня возникли проблемы с генерацией случайной выборки с использованием 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
Спасибо Жан-Филиппу за идею (которая работает).