使用 Python 进行蒙特卡罗模拟：动态构建直方图

Posted 2023-02-16

【中文标题】使用 Python 进行蒙特卡罗模拟：动态构建直方图

【英文标题】：Monte Carlo Simulation with Python: building a histogram on the fly

【发布时间】：2013-08-08 03:20:47

【问题描述】：

我有一个关于使用 Python 动态构建直方图的概念性问题。我正在尝试找出是否有好的算法或现有的包。

我编写了一个函数，它运行蒙特卡罗模拟，被调用 1,000,000,000 次，并在每次运行结束时返回一个 64 位浮点数。以下是上述功能：

def MonteCarlo(df,head,span):
    # Pick initial truck
    rnd_truck = np.random.randint(0,len(df))
    full_length = df['length'][rnd_truck]
    full_weight = df['gvw'][rnd_truck]

    # Loop using other random trucks until the bridge is full
    while True:
        rnd_truck = np.random.randint(0,len(df))
        full_length += head + df['length'][rnd_truck]
        if full_length > span:
            break
        else:
            full_weight += df['gvw'][rnd_truck]

    # Return average weight per feet on the bridge
    return(full_weight/span)

df 是一个 Pandas 数据框对象，其列标记为 'length' 和 'gvw'，分别是卡车长度和重量。 head 是两辆连续卡车之间的距离，span 是桥梁长度。只要卡车列车的总长度小于桥的长度，该功能就会将卡车随机放置在桥上。最后，计算每英尺桥上存在的卡车的平均重量（桥上存在的总重量除以桥长）。

因此，我想构建一个表格直方图，显示返回值的分布，以后可以绘制。我有一些想法：

继续在 numpy 向量中收集返回值，然后在完成 MonteCarlo 分析后使用现有的直方图函数。这是不可行的，因为如果我的计算是正确的，我只需要 7.5 GB 的内存用于该向量（1,000,000,000 64 位浮点数 ~ 7.5 GB）

初始化具有给定范围和 bin 数量的 numpy 数组。在每次运行结束时将匹配箱中的项目数增加一。问题是，我不知道我会得到的值的范围。设置具有范围和适当 bin 大小的直方图是未知的。我还必须弄清楚如何将值分配给正确的 bin，但我认为这是可行的。

以某种方式即时执行。每次函数返回一个数字时修改范围和 bin 大小。我认为这太棘手了，无法从头开始编写。

好吧，我敢打赌可能有更好的方法来处理这个问题。欢迎任何想法！

第二点，我测试了运行上述函数 1,000,000,000 次，只是为了得到计算出的最大值（代码 sn-p 如下）。当span = 200 时，这需要大约一个小时。如果我运行更长的跨度，计算时间会增加（while 循环运行的时间更长以用卡车填充桥梁）。你认为有没有办法优化这个？

max_w = 0
i = 1
    while i < 1000000000:
        if max_w < MonteCarlo(df_basic, 15., 200.):
            max_w = MonteCarlo(df_basic, 15., 200.)
    i += 1
print max_w

谢谢！

【问题讨论】：

给一个 bin 赋值是简单的二分查找。但是，您不能即时更改范围，这意味着您必须提前知道它或存储所有内容。或者至少，做一些假设：例如，将数据聚集在给定大小的小 bin 中（因此您不需要存储太多数据），并在数据“溢出”它们时扩展 bin 列表。

@arbautjc 感谢您的回答。我在最后编辑了与性能问题相关的帖子，但是与我的直方图问题相比，它的优先级较低。我有点希望有一个能够做到这一点的科学软件包。

我给你一个快速而肮脏的实现，使用哈希表而不是排序列表（简单得多）。

【参考方案1】：

这是一个可能的解决方案，具有固定的 bin 大小，并且 bin 的形式为 [k * size, (k + 1) * size[。函数 finalizebins 返回两个列表：一个带有 bin 计数 (a)，另一个 (b) 带有 bin 下限（上限通过添加 binsize 推导出）。

import math, random

def updatebins(bins, binsize, x):
    i = math.floor(x / binsize)
    if i in bins:
        bins[i] += 1
    else:
        bins[i] = 1

def finalizebins(bins, binsize):
    imin = min(bins.keys())
    imax = max(bins.keys())
    a = [0] * (imax - imin + 1)
    b = [binsize * k for k in range(imin, imax + 1)]
    for i in range(imin, imax + 1):
        if i in bins:
            a[i - imin] = bins[i]
    return a, b

# A test with a mixture of gaussian distributions

def check(n):
    bins = 
    binsize = 5.0
    for i in range(n):
        if random.random() > 0.5:
            x = random.gauss(100, 50)
        else:
            x = random.gauss(-200, 150)
        updatebins(bins, binsize, x)
    return finalizebins(bins, binsize)

a, b = check(10000)

# This must be 10000
sum(a)

# Plot the data
from matplotlib.pyplot import *
bar(b,a)
show()