优化用于仿真的Python/Numpy代码
我在做一个模拟项目。为了简单起见:我们的目标是模拟事件的频率，然后模拟每个模拟事件的严重程度，并将这些事件相加一段时间。我试图将其作为一个大小的numpy数组(句点，2)。每一行代表一个时间段。第一个条目包含频率。第二项是对严重程度的总和。我还允许对个人和综合严重程度设定上限。目前，我使用泊松作为频率，帕累托用于严重程度。
我有两种工作方法，但我相信还有更快的方法。目前，如果我模拟10k周期，我的代码需要1.3秒，如果我模拟100 k周期，则需要13秒。这可能取决于我的分布和参数选择，因为更高的频率将意味着更严重的模拟。我最大的努力似乎是在一段时间内对严重程度进行求和/将该总和重新分配到正确的频率行。
这个项目的最终目标是允许多达一百万的模拟多达十个场景，所以更快的速度将是极好的。有小费吗？
尝试1:遍历行
import numpy as np
def simulation(simNum,lambdaPoisson, alphaPareto, claimLimit=np.inf, annualLimit=np.inf):
results = np.empty([simNum, 2])
results[:,0] = np.random.default_rng().poisson(lambdaPoisson, simNum)
for row in results:
row[1] = np.random.default_rng().pareto(alphaPareto, int(row[0])).clip(max=claimLimit,min=None).sum().clip(max=annualLimit,min=None)
print(results)
return results
results = simulation(simNum = 100000, lambdaPoisson = 2, alphaPareto=.2, claimLimit=25, annualLimit = 100)尝试2:映射函数
def mapPareto(a, b, c, d):
out = np.random.default_rng().pareto(b, a).clip(max=c,min=None)
return out.sum().clip(max=d,min=None)
def simulation(simNum,lambdaPoisson, alphaPareto, claimLimit=np.inf, annualLimit=np.inf):
results = np.ones([simNum, 2],dtype=int)
results[:,0] = np.random.default_rng().poisson(lambdaPoisson, simNum)
results[:,1] = list(map(mapPareto, results[:,0], np.full(simNum,alphaPareto),np.full(simNum,claimLimit),np.full(simNum,annualLimit)))
print(results)
return results
results = simulation(simNum = 100000, lambdaPoisson = 2, alphaPareto=.2, claimLimit=25, annualLimit = 100)回答 2
Stack Overflow用户
发布于 2022-07-27 11:51:15
想法1:重用RNG
从重用rng对象开始。而不是
for row in results:
row[1] = np.random.default_rng().pareto(alphaPareto, int(row[0])).clip(max=claimLimit,min=None).sum().clip(max=annualLimit,min=None)使用
for row in results:
row[1] = rng.pareto(alphaPareto, int(row[0])).clip(max=claimLimit,min=None).sum().clip(max=annualLimit,min=None)前(100 K)：9.19秒
后(100 K)：3.09秒
想法2:使用多个核心
而不是
results = simulation(simNum = 1_000_000, lambdaPoisson = 2, alphaPareto=.2, claimLimit=25, annualLimit = 100)使用
def f(_):
return simulation(simNum = 100_000, lambdaPoisson = 2, alphaPareto=.2, claimLimit=25, annualLimit = 100)
pool = multiprocessing.Pool(8)
results = pool.map(f, range(10))前(1M)：30.4秒
后(1M)：9.1秒
在上述两个测试中都使用了Idea 1。
想法3:使用Numba
@DominikStań捷克的回答彻底解释了这一点。
想法4:用矢量化代码替换for
使用一些魔法(我会解释)，您的代码可以重写如下：
import numpy as np
def simulation(simNum,lambdaPoisson, alphaPareto, claimLimit=np.inf, annualLimit=np.inf):
rng = np.random.default_rng()
psizes = rng.poisson(lambdaPoisson, simNum)
samples = rng.pareto(alphaPareto, np.sum(psizes) + 1).clip(max=claimLimit,min=None)
slices = np.cumsum(psizes) - psizes
results = np.add.reduceat(samples, slices).clip(max=annualLimit,min=None)
results[psizes == 0] = 0
print(results)
return psizes, results
results = simulation(simNum = 100_000, lambdaPoisson = 2, alphaPareto=.2, claimLimit=25, annualLimit = 100)前(100 K)：9.19秒
后(100 K)：0.202秒
后(10米)：4.8秒
神奇之处:我使用了一个平面数组来创建用于所有模拟(samples)的Pareto示例。然后，我把求和应用到样本的切片上，为不同的模拟找到和。这是使用reduceat完成的。
抓到：
- 使用了一个额外的虚拟示例，以便
reduceat不会在psizes[-1] == 0的情况下抱怨，因为在这种情况下，我们在slice中将有一个len(samples)的值。虚拟样本的效果将在行results[psizs == 0] = 0.
中被覆盖。
results[psizs == 0] = 0是用来覆盖reduceat在片中存在两个相等值时所采取的默认操作(这是从样本中获取值，而不是返回零)。
结论
这两个想法3和4似乎足够快。Idea 4更紧凑，不需要额外的库，但是Idea 3会产生更易读、更灵活的代码。
Stack Overflow用户
发布于 2022-07-27 09:52:36
下面是一个基于numba的解决方案，可以很好地实现1K目标：
import numpy as np
import numba
@numba.njit
def clip(a, minval, maxval):
if a < minval:
return minval
elif a > maxval:
return maxval
else:
return a
@numba.njit
def simulation3(
simNum, lambdaPoisson, alphaPareto, claimLimit=np.inf, annualLimit=np.inf
):
col0 = np.random.poisson(lambdaPoisson, simNum)
col1 = np.empty(simNum, dtype=np.float64)
for i, row in enumerate(col0):
pareto_events_for_row = np.random.pareto(alphaPareto, row)
clipped = [
clip(pareto_event, maxval=claimLimit, minval=0.0)
for pareto_event in pareto_events_for_row
]
row_sum = sum(clipped)
row_sum_clip = clip(row_sum, maxval=annualLimit, minval=0.0)
col1[i] = row_sum_clip
return np.vstack((col0, col1))
%timeit results = simulation3(simNum = 1_000_000, lambdaPoisson = 2, alphaPareto=.2, claimLimit=25, annualLimit = 100)
# 362 ms ± 16.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)我不得不写我自己的clip，因为np.clip不知何故不同意在numba内剪短单个浮标。
我最初以为您可以在一个大的一维数组中生成col0.sum() Pareto随机变量，然后对第二列中的组和和进行一些奇特的索引，但最终编写起来更简单。我完全希望，如果我们看到一个更快的版本，它将是。
警告-- numba和numpy不共享相同的RNG (这是一个持续的问题)，所以即使您使用np.random.seed，也不会看到完全相同的值。这段代码很简单，适合我的大脑，所以我不认为与你写的代码相比有任何错误，但是随着它的复杂性的增加，最好记住，你只能在统计上比较这两种代码。
https://stackoverflow.com/questions/73135338
复制相似问题
腾讯云开发者
