Подтвердить что ты не робот

Моделирование Монте-Карло с помощью Python: создание гистограммы на лету

У меня есть концептуальный вопрос о построении гистограммы "на лету" с помощью 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 ГБ памяти только для этого вектора (1 000 000 000 64-битных плавающих ~ 7,5 ГБ).

  • Инициализировать массив numpy с заданным диапазоном и количеством ящиков. Увеличьте количество элементов в соответствующем бункере по одному в конце каждого прогона. Проблема в том, что я не знаю диапазон значений, которые я бы получил. Установка гистограммы с диапазоном и соответствующим размером буфера неизвестна. Мне также нужно выяснить, как присвоить значения правильным бункерам, но я думаю, что это выполнимо.

  • Сделайте это как-то на лету. Изменяйте диапазоны и размеры бункера каждый раз, когда функция возвращает число. Это было бы слишком сложно писать с нуля, я думаю.

Хорошо, я уверен, что может быть лучший способ справиться с этой проблемой. Любые идеи приветствуются!

Во второй ноте я протестировал выполнение указанной функции на 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

Спасибо!

4b9b3361

Ответ 1

Вот возможное решение с фиксированным размером бина и бункерами формы [k * size, (k + 1) * size [. Функция finalizebins возвращает два списка: один с подсчетом bin (a), а другой (b) с нижними границами бинов (верхняя граница выводится путем добавления 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()

enter image description here