马尔可夫链蒙特卡洛（python，numpy）

Question

我正在做一些物理方面的研究，需要用一种叫做马尔可夫链蒙特卡洛（MCMC）的方法来分析数据。我尝试自己写一个，但总是遇到问题，因为在使用python和numpy时，某些非常非常小的数字会被四舍五入成零。比如，当我需要计算 numpy.exp(-1000) 这样的表达式时，这个表达式是更大数学公式的一部分，所以我不能简单地对它取对数。

我知道python有一些现成的MCMC模块，我也看了一些，但对它们的文档理解起来有点困难，无法应用。有没有人能推荐一个？我现在有一列数据，要把它放入一个概率分布中。这个分布还有另外两个变量，我会对它们进行随机游走，并记录每一步在马尔可夫链中的情况。然后，我需要根据马尔可夫链为这两个变量制作直方图。如果这个问题太模糊，我很抱歉。任何想法或建议都非常感谢！

Answer 1

用PyMC吧，真的很不错。跟着教程一步步来，你很快就能明白怎么建立一个模型了。

http://pymc-devs.github.io/pymc/tutorial.html

Answer 2

如果你的系统支持更高精度的浮点数，尽量使用它们。例如，如果你有 float128：

import numpy as np
print(np.exp(np.float128(-1000)))  # 5.07595889755e-435
print(np.exp(np.float128(-10000)))  #  1.13548386531e-4343

另外，还可以看看 longdouble。具体支持哪些功能，取决于你的操作系统。

你可以把需要高精度的数组转换成这种格式，然后用 Numpy 的函数来处理它们：

# Example array with 3 dimensions
d = np.random.uniform(-10000, -100, 24)
d.shape = (2, 3, 4)

# Cast to a higher precision
D = d.astype(np.float128)
np.exp(D[:,2])  # array([[4.263772e-4326, 4.3465066e-1474, ...