python - 马尔可夫链蒙特卡罗(python，numpy)

Question

我正在做一些物理学研究，为此我需要使用马尔可夫链蒙特卡罗(MCMC)分析一些数据。我尝试自己写一个，但是当 python/numpy 将一个非常小的数字舍入为零时，我不断遇到错误。特别是当我需要执行诸如 numpy.exp(-1000) 之类的操作时。该表达式本身是一个更大的数学方程的一部分，所以我不能只取它的对数。

我知道有可用于 python 的 MCMC 模块，我已经查看了其中一些模块，但在理解它们的文档以应用它们时遇到了困难。有人可以推荐一个吗？我拥有的是插入概率分布的一列数据。这个分布还有另外两个变量，我将对其进行随机游走并记录马尔可夫链中的每个步骤。然后，我需要根据马尔可夫链为这两个变量中的每一个制作一个直方图。如果这个问题太模糊，我深表歉意。非常感谢任何想法或建议，谢谢!

Answer 1

使用更高精度的 float (如果您的系统可用)。例如，如果您有 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, ...