如何从numba中有效地解压Python中的Monte模拟？解出

Question

我正试图有效地创建蒙特卡罗模拟，因为在我的用例中，我需要运行这个模拟70*10^6次。我希望一个更有经验的人，特别是在表演方面，能给我一些关于我可以尝试什么的想法。我有以下投入：

• Demand
  • 每列是一个产品，每一行是一个月
  • ，一些产品在一个确定的月份中有一个需求由三角形分布元组(最小，平均，最大)估计。对于这些值，我将做一个蒙特卡罗模拟1000 times

• Stock

我想要的输出是查找：

available_products(available_products=stock-demand).之和分布的

• 中值(np.median(np.sum(Available_products)，中值接收1000个模拟的
• 之和。

不过，我有一些问题：

• 的速度，我有直觉，有聪明的方法来计算，利用矢量化的函数。然而，我想不出任何一个，所以我尝试了通常的循环。如果您有任何不同的方法可以更快，请告诉我，
• 固定不能将值设置为数组，在我的解决方案中不能使用demand_jindex_demand_not_ = dict_demand_values_simulationsk
  • 解决方案来设置值，我只需要通过demand_jrow，demand_jindex_demand_not_直接访问demand_j位置

下面是@Glauco建议的使用3D数组满足需求的代码：

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

谢谢

Answer 1

内部循环可以使用数组编程+花式索引来删除，这加快了分配给demand_j的速度。另一点是，您可以生成一次demand_j添加一个维度( num_sim_montecarlo)成为3d数组，而在循环中您必须只读取避免在每个循环中创建值的值。