我试图在python中计算泊松分布,如下所示:
p = math.pow(3,idx)
depart = math.exp(-3) * p
depart = depart / math.factorial(idx)
idx 范围从 0
但我得到了OverflowError: long int too large to convert to float
我试图将离开转换为float
但没有结果。
阶乘变大真的很快:
>>> math.factorial(170)
7257415615307998967396728211129263114716991681296451376543577798900561843401706157852350749242617459511490991237838520776666022565442753025328900773207510902400430280058295603966612599658257104398558294257568966313439612262571094946806711205568880457193340212661452800000000000000000000000000000000000000000L
注意L
; 170 的阶乘仍然可以转换为浮点数:
>>> float(math.factorial(170))
7.257415615307999e+306
但下一个阶乘太大:
>>> float(math.factorial(171))
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
OverflowError: long int too large to convert to float
您可以使用该decimal
模块;计算会更慢,但Decimal()
该类可以处理这种大小的阶乘:
>>> from decimal import Decimal
>>> Decimal(math.factorial(171))
Decimal('1241018070217667823424840524103103992616605577501693185388951803611996075221691752992751978120487585576464959501670387052809889858690710767331242032218484364310473577889968548278290754541561964852153468318044293239598173696899657235903947616152278558180061176365108428800000000000000000000000000000000000000000')
您必须在Decimal()
整个过程中使用值:
from decimal import *
with localcontext() as ctx:
ctx.prec = 32 # desired precision
p = ctx.power(3, idx)
depart = ctx.exp(-3) * p
depart /= math.factorial(idx)
当idx
变大时,math.pow
和/或math.factorial
将变得非常大并且无法转换为浮点值(idx=1000
在我的 64 位机器上触发错误)。您不希望使用 math.pow 函数,因为它比内置**
运算符更早溢出,因为它试图通过更早的浮点转换来保持更高的精度。此外,您可以将每个函数调用包装在一个Decimal
对象中以获得更高的精度。
处理非常大的数字时,另一种方法是使用对数刻度。取每个值的对数(或计算每个值的对数版本)并在对结果取幂之前执行所有必需的操作。这允许您的值暂时离开浮动域空间,同时仍然准确计算位于浮动域内的最终答案。
3 ** idx => math.log(3) * idx
math.exp(-3) * p => -3 + math.log(p)
math.factorial(idx) => sum(math.log(ii) for ii in range(1, idx + 1))
...
math.exp(result)
这会一直保留在日志域中,直到最后,所以在遇到溢出问题之前,您的数字会变得非常非常大。
尝试使用十进制库。它声称支持任意精度。
from decimal import Decimal
此外,您不需要使用math.pow
. pow
是内置的。
该scipy
模块可以帮助您。
scipy.misc.factorial是一个阶乘函数,可以使用伽玛函数逼近计算阶乘,并使用浮点返回结果。
import numpy
from scipy.misc import factorial
i = numpy.arange(10)
print(numpy.exp(-3) * 3**i / factorial(i))
给出:
[ 0.04978707 0.14936121 0.22404181 0.22404181 0.16803136 0.10081881
0.05040941 0.02160403 0.00810151 0.0027005 ]
还有一个计算泊松分布的模块。例如:
import numpy
from scipy.stats import poisson
i = numpy.arange(10)
p = poisson(3)
print(p.pmf(i))
给出:
[ 0.04978707 0.14936121 0.22404181 0.22404181 0.16803136 0.10081881
0.05040941 0.02160403 0.00810151 0.0027005 ]