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

Question

我有一个关于使用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是桥长.只要卡车列车的总长度小于桥梁长度,该功能就会随意将卡车放在桥上.最后,计算每英尺桥上存在的卡车的平均重量(桥上存在的总重量除以桥长度).

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

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

2. 初始化具有给定范围和数量的bin的numpy数组.每次运行结束时,将匹配区域中的项目数增加1.问题是,我不知道我会得到的价值范围.设置具有范围和适当的箱尺寸的直方图是未知的.我还必须弄清楚如何将值分配给正确的箱子,但我认为它是可行的.

3. 以某种方式做它在飞行中.每次函数返回一个数字时,修改范围和bin大小.我认为这从头开始写起来太棘手了.

好吧,我打赌可能有更好的方法来处理这个问题.任何想法都会受到欢迎!

在第二个注释中,我测试运行上述函数1,000,000,000次只是为了获得计算的最大值(下面的代码片段).这需要大约一个小时的时间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

谢谢!

Answer 1

这是一个可能的解决方案，具有固定的 bin 大小，以及 [k * size, (k + 1) * size[ 形式的 bin。函数 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()