[閒聊] 替RandomVariate提供提取seed的功能

[AmibaGelos (Amiba Gelos)]

最近在用MMA跑Monte Carlo,覺得MMA不能提取seed實在是很不方便,所以就做了一個可以
提取seed的RandomVariate

Module[{seed, store = {}, temp = 1, Storing},
  RandomMemory[1, opt___] := RandomMemory[temp, opt];
  RandomMemory[in__] := (Storing[in]; RandomVariate[in]);
  Storing[dist_, opt___] :=
   AppendTo[
    store, {Evaluate[If[dist === temp, 1, temp = dist]], opt}];
  SetSeed[InSeed_, InStore_: {}] := (store = {}; temp = 1;
    SeedRandom[seed = InSeed]; (RandomMemory[##];) & @@@ InStore;);
  ExtractSeed[] := {seed, store}; SetSeed[0];];

這個設計在跑custom distribution的時候蠻好用的, 不過AppendTo會是瓶頸

目前只有想到hash, 不過有新distribution還是得重建table

不知道有沒有更漂亮的寫法?

--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 140.112.4.208
※ 文章網址: https://www.ptt.cc/bbs/Mathematica/M.1476160101.A.361.html
→ ginstein: 不懂! seed不就使用者設定的嗎? 不能提取自己設的seed?
推 ginstein: 不懂!seed不就使用者設定的嗎?不能提取自己設的seed?
推 LPH66: 他要的是跑到一半記錄目前 PRNG 的 state 吧
→ LPH66: 不過記得不少演算法的 state 並不只有 seed 所以有點難
→ LPH66: 提供一個統一界面, 例如 MT19937 需要 624 個 32 位元整數
→ LPH66: 但種子只有 32 位元, 所以不可能用單一整數表示狀態
→ LPH66: 然後, 剛剛翻 Help 有看到 BlockRandom 或許也符合你的需要
→ AmibaGelos: 謝謝樓上解釋~MMA預設CA是80*64b,最小是64*64所以不可
→ AmibaGelos: 能用long描述, MMA的seed我覺得實際上是被hash處理過
→ AmibaGelos: 這方面我還在測試,不過想從這邊retrive state也是很難
→ AmibaGelos: BlockRandom只是freeze state而已,應該不能做到同時跑
→ AmibaGelos: 兩個prng,更做不到resume吧@@
→ AmibaGelos: 其實有種最偷懶的做法是每次要存檔的時候就另開一個
→ AmibaGelos: kernel專門跑rng,用mathlink連到另一個MMA也許更clear
→ AmibaGelos: 忽然覺得這應該就是最佳解了@@
→ AmibaGelos: 唯一可惜的是還是得額外存一套distribution,跑MCMC有
→ AmibaGelos: 點浪費ram...
→ ginstein: 謝謝L兄,懂了.另MMAseed代表所有狀態,沒看過找seed指令