※ 引述《[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海