马尔可夫链蒙特卡罗（python，numpy）答案

【问题标题】：Markov Chain Monte Carlo (python, numpy)马尔可夫链蒙特卡罗（python，numpy）
【发布时间】：2014-09-26 08:20:55
【问题描述】：

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

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

【问题讨论】：

您到底是如何计算出某些内容被舍入为 0 的？您是否尝试将 numpy 数字类型的精度提高到 64 位？
我知道它正在向下舍入为 0，因为我收到错误：RuntimeWarning: 除以零在日志中遇到。我正在尝试做类似 np.log(np.exp(-2000)) 的事情。我已经尝试过 np.float128() 方法，但我传递了一个数字数组，而 float 128 一次只能做一个。
你可以只处理概率的日志吗？

标签： python numpy montecarlo markov-chains

【解决方案1】：

如果您的系统上可用，请使用更高精度的浮点数。例如，如果您有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, ...

【讨论】：

感谢您的建议，它几乎解决了我的问题。问题是我不会只将一个数字传递给浮点数。我将一维数组/列传递给它。更多地把它想象成 np.exp(np.float128(d[:,3]))。我需要对整个数组进行计算。有没有办法一次为整个数组做这样的事情？或者我真的需要逐个元素地做吗？谢谢
@user1679198 查看更新后的答案astype 将数组d 转换为更高精度的数组D
您不能使用日志转换并避免exp 调用吗？

【解决方案2】：

使用 PyMC - 很棒。完成教程，您应该很快就会了解如何构建模型。

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

【讨论】：