※ 引述《[email protected] (野孩子)》之銘言： > 這是我為了算negative binomial distribution所寫的程式 > 我是算出P(y=i),i=1...146 > 再由10000個等份中找出他所占的等份 > 簡單來說 > 就是要跑出一份符合分佈的亂數 > 但很不幸的是 > 本來應該E(x)=24,var(x)=600 > 但我程式跑出來的結果卻是E(x)=21,var(x)=510 (約略) > 我一直找不到我錯誤的地方在哪 > 所以只好來請教各位大大 > 希望知道我錯在哪的可以跟我說一下 > 真的感激不盡Orz... > #include<stdio.h> > #include<stdlib.h> > #include<time.h> > long int factorial(int p); > double power(double,int); > long int factorial(int p) > { > int count; > long int result=1; ^^^^^^^^^^^^^^^^^^ 對 c 不熟, 不確定這樣寫會不會造成問題. 建議改用指令 result=1; 以保證每次呼叫時, result 都有正確起始值 1. > for(count=1;count<=p;count++) > { > result=result*count; > } > return result; > } 檢查程式, 尤其是別人寫的長程式, 是很耗時耗精神的工 作, 所以以下我不再仔細檢查. 但有幾點可能要注意: (1) 電腦儲存數字的精確度有限, 以連乘、連加計算時, 要注意是否會有很大的累積誤差。 (2) 計算二項係數涉及的數字很大, 可能造成 overflow, 可考慮用對數。 (3) 計算一個小數的乘冪, 若乘冪太大, 要注意是否會造 成 underflow. (4) 負二項 r.v. 可以由幾何 r.v. 建構, 而幾何 r.v. 較容易做(機率計算容易, 且不需累加). 因此, 若應 用上不需太高的 "成功數", 可考慮此法。 (5) 在成功機率 p 不是很小情況, 可模擬 Bernoulli 實 驗以產生所要的負二項變量, 或許更快? (6) 如果 "寫產生亂數的程式" 不是望要目的, 可網路上 找找, 或許有現成副程式庫可用.(商業軟體IMSL有) -- ◢ 川◣ │││││ 您在找統計版嗎? 竭誠邀請您加入 Statistics! ▃▅▃▅▆ ◣││││ 無名小站 telnet://wretch.twbbs.org (cat▍_/ ▲ 、 ││ 成大計中站 telnet://bbs.ncku.edu.tw ▊ ▊Moon◤▍▍▄▂ │ 交大次世代 telnet://bs2.twbbs.org ▃─ _▍_ ◣▌▎▃▅ 盈月與繁星 telnet://ms.twbbs.org ▲ ◤ ￣ ◢▂▃ ＊Mooncat~ ★未經本人同意請勿轉載; 回覆請勿全文引用! -- 夫兵者不祥之器物或惡之故有道者不處君子居則貴左用兵則貴右兵者不祥之器非君子 之器不得已而用之恬淡為上勝而不美而美之者是樂殺人夫樂殺人者則不可得志於天下 矣吉事尚左凶事尚右偏將軍居左上將軍居右言以喪禮處之殺人之眾以哀悲泣之戰勝以 喪禮處之道常無名樸雖小天下莫能臣侯王若能守之萬物將自賓天地相合以降甘露民莫 之令而自均始制有名名亦既有夫亦將知止知止可以不殆譬道之在天 163.15.188.87海