本文摘於 Orchis(梅之雅在雪前香) Orchis.bbs@cis.nctu.edu.tw ddj(聖彼得逆鉤參) ddj.bbs@bbs.iie.ncku.edu.tw Devil(璉璉) Devil.bbs@bbs.iie.ncku.edu.tw yhliu(劉應興) yhliu.bbs@bbs.iie.ncku.edu.tw GGQ(一杯酒情緒萬種.) GGQ.bbs@bbs.iie.ncku.edu.tw ---------------------------------------------------------------------------- 發信人: Orchis.bbs@cis.nctu.edu.tw (梅之雅在雪前香), 信區: fortran 標 題: 亂數產生器.... 請問有一種方法其所產生的亂數是呈現均勻分佈 的嗎??如: 我們要50個數...而其會均勻分佈在 0與1之間... 謝謝.. ---------------------------------------------------------------------------- 發信人: ddj@bar (聖彼得逆鉤參), 信區: fortran 標 題: Re: 亂數產生器.... 若以MS-Fortran 5.1為例: integer*2 ihr,imin,isec,ipsec ......... call gettim(ihr,imin,isec,ipsec) --- 叫用系統時間 call seed(ipsec) ----以百分之一秒為亂數產生器起始種子 do i=1,50 call random(F) ----產生0~1之均勻亂數(F值) ....... end do ps:但有一個問題喔!每次執行時產生的第一個亂數值都很小,每次ㄛ! 是不是上述ipsec是指程式執行後至該行所花時間中百分秒數,所 以每次執行時都差不多...??? 若是,怎麼改善?若不是,why?? *若在工作站跑,寫法又不一樣,可參考使用手冊.... 發信人: Devil@bar (璉璉), 信區: fortran 補充說明: 1.目前電腦上沒有所謂的真亂數及均勻分佈, 但可以非常近似, 本版網友有統計 所的人材, 相信他在這部份可以給各位網友更詳細的解釋, 有興趣者可自行詢問 !!當然, 他想主動提出的話, 我一定會收錄在精華區的!! 2.目前的亂數通常都採用一個很大的質數再取餘數, 所以取的亂數量很大時一定 會造成循環, 這時就可能要考慮自行建造亂數產生器. 3.由於上述程式碼的亂數種子介於 0-99 , 可以考慮較大的亂數種子. 如 call seed(ipsec*100+isec) 也許可得到較佳的結果. ---------------------------------------------------------------------------- 發信人: yhliu@bar (劉應興), 信區: fortran 標 題: Re: 亂數產生器.... 關於亂數產生器, 我所知也很有限。因為 Devil 版主提到有 統計所的人 (不一定是指我 ^_^), 所以我也來野人獻曝。 Devil 說得很對, 目前常用的亂數, 多是以「乘積同餘法」 (multiplicative-congruential method) 做的; 一個稍為好 一點的方法是 congruential method, 即 x(i) = (a*x(i-1)+c) (mod m) 但多加一個常數, 增加運算的時間。雖然就產生少數亂數而 言, 那是不重要的, 但對經常需要極大量亂數的計算而言, 事實上產生亂數所需時間是很重要的。有 IMSL 的 source 的網友不妨 check 一下, IMSL 是否有加上一個常數﹖為了 節省計算時間, 上述 m 也常用 2 的乘冪, 如古老的 IBM SSP 副程式庫, 就是採用 2^{16} 或 2^{32} (視 2 bytes 整數或 4 bytes 整數而定), 其原理是根本不用再做同餘運 算 (overflow 自動去掉了)。Kennedy and Gentle (1980) 的 "Statistical Computing" p.145 有抄錄一個用 IBM S/360 組合語言寫的副程式 (副程式本體只有五個運算指令)。 事實上「同餘法」的均勻性並不好, 這方面的討論可參考上 述 Kennedy and Gentle (1980) 的書, 或 Knuth (1969), "The Art of Computer Programming", Vol.2。這裡所講的 「均勻性」不是只看一維的, 而是看連續幾個亂數的分布; 如二維的情形, 前引書有用簡單的例子顯示這樣的問題。當 然除數較大情形可能較不嚴重, 不過, 基本的問題還是存在 的。Tausworthe (1965) 建議了一種方法, 稱為 feedback shift register (FSR) 的亂數產生器, Lewis and Payne (1973) 做了一些修改, 稱為 GFSR (generalized FSR), 此法產生 的亂數, 其均勻性似乎較好。在 Kennedy and Gentle 的書 pp.160-161 有完整的副程式。我曾試過, 看二維的分布, 看 起來相當好。 事實上, 除非加上一些程式外的資料, 否則任何亂數產生器 一定是循環的, 只差循環週期的長短。原因是亂數產生器是 一個固定的運算程序, 而其前一個輸出是下一個運算的輸入, 但計算機的數值儲存狀態是有限的, 因此必然會產生與前面 相同的結果, 則開始一個新的循環週期。這有點像有理數表 示為小數必然是循環的一樣。要延長週期, 可採用兩個以上 的亂數產生器。有許多方式可組合應用兩個產生器, 例如用 第一個產生器產生 k 個亂數, 用第二個亂數產生器決定選用 其中一個為輸出之亂數; 或者, 可以用第一個亂數產生器產 生一個 seed, 並產生一個整數 k, 而用第二個產生器. 以 seed 為種子, 產生 k 個亂數。 好的亂數產生器其輸出的均勻性應該不會受 seed 很大的影 響 (以乘積同餘法來說, 重要的是建立在程式中的乘數和除 數)。但取一個比較大的 seed 是合理的, 因為這樣我們的 第一個 output 才有點「隨機」的味道。或者, 可把前幾個 output 丟棄, 也就是主程式中先做一個 loop, 呼叫幾次亂 數產生器, 然後以其輸出為 seed 開始做模擬。 產生亂數的方法, 目前還有很多人在做研究。我不是專長這 方面的, 以上所述的主要是參考上述兩本書----事實上主要 是第一本 (因 Knuth 的書沒帶來)。沒有甚麼學問, 卻寫得 很嚕嗦, 請看官們原諒! ^_^ ---------------------------------------------------------------------------- 發信人: GGQ@bar (一杯酒情緒萬種.), 信區: fortran 發信站: 成大資訊所_BBS ( Sep 20 19:48:24 1998) 我再次去查看IMSL的使用手冊. 這次終於看懂關於RNSET(),和RNGET()的部份了. 也就是說,如果你第一次是用 call rnset(0) 的話,那麼它第一次就會以system clock來算,作為種子(iseed) 之後,你只要再call rnget(iseed)的話, 它就會以第一次的亂數中再去獲的另一個新的種子(iseed) 此時它會以iseed這一個變數作為OUTPUT喔. 那你就再去 call rnset(iseed) 也就是說,你又告訴了電腦說,這一次你要以iseed為種子了喔. 如此,如果你想利用這兩個函式庫RNSET(),和RNGET(),那麼就並不算難了. 不過,要注意的一點是,RNGET(),跟RNSET(), 只適用在,IMSL的random number generator.才會有作用到. 不然,你會發現,可能種子沒再便喔. 舉個例來說,如果我想產生十次的種子,讓它都不相同的話. 順便用IMSL裡面的其中一項RNUNF去產生十個uniform(0,1)的數字,順便 把它的種子也印出來看看. program GGQ integer i,iseed real R(10),RNUNF external rnset,rnget,rnunf iseed=0 <-我第一次設iseed=0 do i=1,10 call rnset(iseed) <--此時告訴電腦,我以iseed為INPUT當作種子個數,當然第一次 iseed=0,也就是會以system clock來產生種子,當然之後 就不會了,iseed會一直去改變. R(i)=RNUNF() <--這裡用到IMSL裡面的RNUNF()去產生服從UNiform(0,1)的亂數 write(*,*) iseed <--把種子印出來看看. call rnget(iseed) <--此時告訴電腦我還要去,產生另一個種子,名為iseed的變項. 是以iseed為OUTPUT end do write(*,*) 'Uniform(0,1)'s 10 random number',R end _________________________________________________________________ 真有趣喔,好像比那個我原本用的方法好. call gettim(ihr,imin,isec,ipsec) iseed=100000*ipsec ________________________________ 作為大家參考,我也是走一步學一步,大家參考切搓罷.