Генерация одинаковых случайных чисел в R и Julia

Я хотел бы генерировать одинаковые случайные числа в R и Julia. Оба языка по умолчанию используют библиотеку Mersenne-Twister, однако в Julia 1.0.0:

julia> using Random
julia> Random.seed!(3)
julia> rand()
0.8116984049958615

Производит 0.811..., а в R:

set.seed(3)
runif(1)

производит 0.168.

Любые идеи?

Связанные вопросы SO здесь и здесь.

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


r random julia mersenne-twister
person Colin T Bowers    schedule 07.04.2015    source источник
comment
Грубым способом было бы создать все реплики начальной загрузки (или, возможно, только индексы) заранее и сохранить их в файле, который могли бы использовать обе программы.   -  person joran    schedule 07.04.2015
comment
Это не ответ, но я предполагаю, что способ перевода семени в начальное состояние для библиотеки MT отличается. Я предполагаю, что ответы могут и должны быть найдены в исходном коде (ура для открытого исходного кода).   -  person IainDunning    schedule 07.04.2015
comment
@joran Согласен, и это то, что я могу в конечном итоге сделать. Однако для этого есть немного работы (по крайней мере, для меня - я относительный новичок в R), поскольку это подразумевает изменение исходного кода R и Julia для поиска случайных чисел в файле.   -  person Colin T Bowers    schedule 07.04.2015
comment
@IainDunning Звучит разумно для меня. Я подумал, что сначала спрошу здесь на всякий случай, если кто-то сможет ответить за 5 минут, что может занять у меня целый день :-)   -  person Colin T Bowers    schedule 07.04.2015
comment
Использование RCall не помогает?   -  person Khashaa    schedule 07.04.2015
comment
@Khashaa RCall, безусловно, помогает мне передавать данные между ними (например, если я генерирую случайный вектор чисел, который я хочу использовать в качестве источника случайности на обоих языках), но, как и в случае с предложением Джорана, это по-прежнему подразумевает, что мне нужно будет отредактировать интересующий исходный код R, чтобы указать на этот случайный вектор. По общему признанию, это именно то, что я, вероятно, в конечном итоге сделаю :-)   -  person Colin T Bowers    schedule 07.04.2015
comment
К вашему сведению, srand(3); rand() на моей 32-битной платформе Linux Julia производит 0.8116984049958615.   -  person rickhg12hs    schedule 08.04.2015
comment
@ rickhg12hs У меня 64-разрядная версия Ubuntu 14.04. Вероятно, следовало упомянуть об этом в вопросе. Но я предполагаю, что эта разница просто подчеркивает смысл ответа Дирка - здесь нет волшебного короткого пути...   -  person Colin T Bowers    schedule 10.04.2015

Ответы (3)


arrow_upward
7
arrow_downward

Это старая проблема.

Пол Гилберт обращался к той же проблеме в конце 1990-х (!!), когда пытался утверждать, что симуляции в R (тогда еще новичке) давали тот же результат, что и в S-Plus (тогда господствующем).

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

person Dirk Eddelbuettel    schedule 07.04.2015
comment
Другими словами, никакого волшебного ярлыка :-) Спасибо за ответ, я воздержусь от постановки галочки на день или два, если кто-то сможет придумать что-то, что сработает в моем конкретном случае. - person Colin T Bowers; 07.04.2015

arrow_upward
5
arrow_downward

Следуя предложению RCall, сделанному @Khashaa, ясно, что вы можете установить начальное число и получить случайные числа от R.

julia> using RCall

julia> RCall.reval("set.seed(3)")
RCall.NilSxp(16777344,Ptr{Void} @0x0a4b6330)

julia> a = zeros(Float64,20);

julia> unsafe_copy!(pointer(a), RCall.reval("runif(20)").pv, 20)
Ptr{Float64} @0x972f4860

julia> map(x -> @printf("%20.15f\n", x), a);
   0.168041526339948
   0.807516399072483
   0.384942351374775
   0.327734317164868
   0.602100674761459
   0.604394054040313
   0.124633444240317
   0.294600924244151
   0.577609919011593
   0.630979274399579
   0.512015897547826
   0.505023914156482
   0.534035353455693
   0.557249435689300
   0.867919487645850
   0.829708693316206
   0.111449153395370
   0.703688358888030
   0.897488264366984
   0.279732553754002

и от R:

> options(digits=15)
> set.seed(3)
> runif(20)
 [1] 0.168041526339948 0.807516399072483 0.384942351374775 0.327734317164868
 [5] 0.602100674761459 0.604394054040313 0.124633444240317 0.294600924244151
 [9] 0.577609919011593 0.630979274399579 0.512015897547826 0.505023914156482
[13] 0.534035353455693 0.557249435689300 0.867919487645850 0.829708693316206
[17] 0.111449153395370 0.703688358888030 0.897488264366984 0.279732553754002

** ИЗМЕНИТЬ **

Согласно предложению @ColinTBowers, вот более простой/чистый способ доступа к R случайным числам из Julia.

julia> using RCall

julia> reval("set.seed(3)");

julia> a = rcopy("runif(20)");

julia> map(x -> @printf("%20.15f\n", x), a);
   0.168041526339948
   0.807516399072483
   0.384942351374775
   0.327734317164868
   0.602100674761459
   0.604394054040313
   0.124633444240317
   0.294600924244151
   0.577609919011593
   0.630979274399579
   0.512015897547826
   0.505023914156482
   0.534035353455693
   0.557249435689300
   0.867919487645850
   0.829708693316206
   0.111449153395370
   0.703688358888030
   0.897488264366984
   0.279732553754002
person rickhg12hs    schedule 12.04.2015
comment
Ницца. Таким образом, есть ярлык через RCall, который, вероятно, можно было бы обернуть. Но это лишь подчеркивает главное: если вам нужен один и тот же поток ГСЧ, вы действительно хотите, чтобы выполнялся один и тот же код. Я мог бы начать с простого генератора (скажем, KISS Марсальи) и просто закодировать их с обеих сторон. - person Dirk Eddelbuettel; 12.04.2015
comment
@DirkEddelbuettel, я искал Mersenne-Twisters с открытым исходным кодом и нашел примеры по адресу Веб-сайт Макото Мацумото (множество версий исходного кода для загрузки и исходный документ, содержащий исходный код), исходный код R и GSL. Все они немного разные. К счастью, C-интерфейс Джулии работает хорошо, а R предоставляет общую библиотеку и т. д. - person rickhg12hs; 12.04.2015
comment
МТ слишком сложный. Используйте что-то вроде простых конгруэнтных линейных генераторов (для создания униформы; например, что-то вроде KISS) в, например, Ziggurat (который создает нормали). Зиккурат используется в Джулии, а у меня есть пакет RcppZiggurat для R. - person Dirk Eddelbuettel; 12.04.2015
comment
Кроме того, глобальный генератор случайных чисел Джулии по умолчанию является последней библиотекой dSFMT. Base.Random.globalRNG() |> typeof возвращает Base.Random.MersenneTwister. R's Mersenne Twister лишь немного изменен по сравнению с оригинальной бумагой. Я думаю, что раздача R тоже битый. - person rickhg12hs; 12.04.2015
comment
Спасибо за ответ. Это интересно, хотя я ошибаюсь, думая, что это немного сложнее, чем должно быть? Не могли бы вы просто использовать (от Джулии) reval("set.seed(3)"), а затем x = rcopy("runif(20)"), чтобы импортировать числа в локальную переменную x в Джулии? Возможно, я что-то упускаю? - person Colin T Bowers; 13.04.2015
comment
@ColinTBowers Вполне возможно. Я RCall нуб и только что сообщил о первом, что я попробовал, и это сработало. Если есть более простой и/или более надежный способ получить R случайных чисел в Julia, я только за! - person rickhg12hs; 13.04.2015

arrow_upward
2
arrow_downward

Видеть:

?set.seed

«Мерсенн-Твистер»: из Мацумото и Нисимуры (1998). Скрученная GFSR с периодом 2^19937 - 1 и равнораспределением в 623 последовательных измерениях (за весь период). «Семя» — это 624-мерный набор 32-битных целых чисел плюс текущая позиция в этом наборе.

И вы можете увидеть, можете ли вы ссылаться на один и тот же код C из обоих языков. Если вы хотите увидеть список/вектор, введите:

.Random.seed
person IRTFM    schedule 07.04.2015
comment
Да, но если бы это было так просто, вы бы также получили те же результаты с GSL Mersenne Twister или другим. Обычно мешают небольшие локальные различия в создании наборов, манипулировании и т. д. Я бы просто написал простую процедуру... - person Dirk Eddelbuettel; 07.04.2015
comment
Если кто-то задается вопросом, кому из нас верить,... верьте Дирку. Вероятно, это сэкономит вам много времени. - person IRTFM; 07.04.2015