Блог

Я пытаюсь эффективно создать симуляцию Монте-Карло, потому что в моем случае мне нужно было бы запустить это моделирование 70 * 10 ^ 6 раз. Я надеялся, что кто-то более опытный, особенно в производительности, может дать мне несколько идей о том, что я мог бы попробовать. У меня есть следующие входные данные:

• Спрос
  • Каждый столбец — это продукт, каждая строка — месяц
  • Спрос на некоторые продукты в определенный месяц оценивается с помощью треугольного кортежа распределения (min, mean, max). Для этих значений я проведу моделирование методом Монте-Карло 1000 раз
• Запас

Мой желаемый результат — найти:

• Медиана распределения суммы доступных продуктов (np.median(np.sum(available_products))), медиана получает 1000 симуляций суммы доступных продуктов (available_products = спрос на акции).

Однако у меня возникли некоторые проблемы:

• Скорость, у меня есть интуиция, что существуют умные способы вычисления векторизованных функций с использованием векторизованных функций. Однако я не мог придумать ни одного, поэтому я попробовал обычные циклы. Пожалуйста, если у вас есть какие-либо подсказки о каком-либо другом подходе, который может быть быстрее, дайте мне знать.
• ИСПРАВЛЕНО Не удается установить значения в массив, в моем решении я не могу установить значения, используя demand_j[index_demand_not_0][k] = dict_demand_values_simulations[k][j]
  • РЕШЕНИЕ, мне просто нужно получить прямой доступ к местоположению demand_j с помощью demand_j[строка, столбец].

Вот код, использующий 3D-массив для спроса, предложенный @Glauco:

 import numpy as np
from numba import jit


@jit(nopython=True, nogil=True, fastmath=True)
def calc_triangular_dist(demand_distribution, num_monte):
    # Calculates triangular distributions
    return np.random.triangular(demand_distribution[0], demand_distribution[1], demand_distribution[2], size=num_monte)


def demand3d():
    # Goal find distribution_of_median_of_sum_available_products(np.median(np.sum(available_products)), the median from the 1000 Monte Carlo Simulations ): available_products=stock-demand (Each demand is generated by a Monte Carlo simulation 1000 times, therefore I will have 1000 demand arrays and consequently I will have a distribution of 1000 values of available products)
    # Input
    demand_triangular = np.array(
        [
            [0.0, 0.0, 0.0, 0.0],
            [0.0, 0.0, 0.0, (4.5, 5.5, 8.25)],
            [(2.1, 3.1, 4.65), 0.0, 0.0, (4.5, 5.5, 8.25)],
        ]
    )  # Each column represents a product, each row a month. Tuples are for triangular distribution (min,mean,max)
    stock = np.array(
        [[30, 30, 30, 22], [30, 30, 30, 22], [30, 30, 30, 22]]
    )  # Stock of available products, Each column represents a product, each row a month.
    num_sim_monte_carlo = 1000

    # Problem 1) How to unpack effectively each array of demand from simulation? Given that in my real case I would have 70 tuples to perform the Monte Carlo simulation?

    row, col = demand_triangular.shape
    index_demand_not_0 = np.where(
        demand_triangular != 0
    )  # Index of values that are not zeros,therefore my tuples for triangular distribution

    demand_j = np.zeros(shape=(row, col,num_sim_monte_carlo), dtype=float)

    triangular_len = len(demand_triangular[index_demand_not_0])  # Length of rows to calculate triangular
    for k in range(0, triangular_len):  # loop per values to simulate
        demand_j[index_demand_not_0[0][k], index_demand_not_0[1][k]] = calc_triangular_dist(
            demand_triangular[index_demand_not_0][k], num_sim_monte_carlo
        )

    sums_available_simulations = np.zeros(
        shape=num_sim_monte_carlo
    )  # Stores each 1000 different sums of available, generated by unpacking the dict_demand_velues_simulations

    for j in range(0, num_sim_monte_carlo):  # loop per number of monte carlo simulations
        available = stock - demand_j[:,:,j]
        available[available < 0] = 0  # Fixes with values are negative
        sums_available_simulations[j] = np.sum(available)  # Stores available for each simulation
    print("Median of distribution of available is: ", np.median(sums_available_simulations))

if __name__ == "__main__":
    demand3d()
 

Результаты предложений показывают гораздо лучшую производительность при использовании 3D-массива :), теперь, когда у меня есть только массивы, я могу попытаться улучшить дальнейшее использование numba.

 Baseline  0.4067141000000001
1) Monte Carlo per loop  0.035586100000000176
2) Demand 3D  0.017964299999999822
 

Спасибо

внутренний цикл можно удалить с помощью программирования массива причудливого индексирования, это ускоряет присвоение demand_j . Другой момент заключается в том, что вы можете сгенерировать один раз demand_j, добавив измерение (num_sim_montecarlo), чтобы оно стало 3D-массивом, и в цикле вы должны читать только значения, избегая создания значений в каждом цикле.