numpy Python -加速for循环

3qpi33ja  于 2023-11-18  发布在  Python
关注(0)|答案(1)|浏览(133)

因此,我目前正在编写一个脚本,它执行以下步骤:
1.根据泊松分布进行模拟。假设此样本根据泊松分布得到5
1.根据泊松分布的值,模拟对数正态分布的5
1.重复上述过程100 k+次
1.对这些值应用一些逻辑以获得一些统计信息
我认为最快的方法是通过一个大型矩阵,其中包含从对数正态分布模拟的每个行和列的每个值。有了这个矩阵,我就可以应用我需要的度量。
这是我到目前为止所做的,这似乎是可行的-然而,我是python的初学者,我认为应该有一种方法可以做到这一点“矢量化”,而不必诉诸于for循环。
任何速度上的改进都将受到极大的赞赏,因为我理想地想运行1 m模拟:

import pandas as pd
import numpy as np
import os
import scipy.stats as sp

n_sims = 100000
MU = 3
poisson = sp.poisson(mu=MU)

#first want to sample 10000 numbers from a poisson ## this is step 1
no_clms_arr = poisson.rvs(size=n_sims)

#now for each element sampled sample a lognormal and store the results back to a matrix
result = np.zeros((n_sims,no_clms_arr.max()))

for idx,i in enumerate(no_clms_arr):
    if i == 0:
        continue # no need to do anything here
    sev_sample = sp.lognorm(50,25).rvs(i) ## this is step 2
    sev_sample = np.pad(sev_sample,(0,no_clms_arr.max()-i),'constant',constant_values=(0,0))  ## padding out so I can append back to the result matrix

    result[idx,:] = sev_sample

#now calculate some metrics ## step 4, example metric is below
lim = 5e6
xs = 2e6
np.clip(a = result-xs,a_min=0,a_max=lim).sum(axis=1).mean()

字符串

uurity8g

uurity8g1#

你可以把它完全矢量化,因为result数组是不必要的,所以你只需要生成你的随机数,然后进行裁剪和求和:

vals = sp.lognorm(50, 25).rvs(no_clms_arr.sum())
(vals - xs).clip(0, lim).sum() / n_sims

字符串
你可以通过一个for循环来加速你的方法(假设有必要使数组result),方法是在一个单独的调用中在for循环之外生成所有的随机数,然后用一个for循环将这些值分配给result中的相应点,即用以下内容替换你的for循环:

vals = sp.lognorm(50, 25).rvs(no_clms_arr.sum())

start = 0

for i, j in enumerate(no_clms_arr):
    if j:
        result[i, :j] = vals[start: start + j]
        start += j


您还可以将.sum(axis=1).mean()调用替换为.sum()/n_sims,以获得额外的速度提升。
时间:

def op(n_sims = 10000, mu = 3, lim=5e6, xs=2e6):
    poisson = sp.poisson(mu=mu)
    no_clms_arr = poisson.rvs(size=n_sims)
    result = np.zeros((n_sims,no_clms_arr.max()))
    for idx,i in enumerate(no_clms_arr):
        if i == 0:
            continue # no need to do anything here
        sev_sample = sp.lognorm(50,25).rvs(i) ## this is step 2
        sev_sample = np.pad(sev_sample,(0,no_clms_arr.max()-i),'constant',constant_values=(0,0))  ## padding out so I can append back to the result matrix
        result[idx,:] = sev_sample
    return np.clip(a = result-xs,a_min=0,a_max=lim).sum(axis=1).mean()

def nin17_vectorized(n_sims = 10000, mu=3, lim=5e6, xs=2e6):
    poisson = sp.poisson(mu=mu)
    no_clms_arr = poisson.rvs(size=n_sims)
    vals = sp.lognorm(50, 25).rvs(no_clms_arr.sum())
    return (vals - xs).clip(0, lim).sum() / n_sims

def nin17_for(n_sims = 10000, mu=3, lim=5e6, xs=2e6):
    poisson = sp.poisson(mu=mu)
    no_clms_arr = poisson.rvs(size=n_sims)
    result = np.zeros((n_sims,no_clms_arr.max()))
    vals = sp.lognorm(50, 25).rvs(no_clms_arr.sum())
    start = 0
    for i, j in enumerate(no_clms_arr):
        if j:
            result[i, :j] = vals[start: start + j]
            start += j
    return np.clip(a = result-xs,a_min=0,a_max=lim).sum() / n_sims

state = np.random.get_state()

np.random.set_state(state)
print(op())
np.random.set_state(state)
print(nin17_for())
np.random.set_state(state)
print(nin17_vectorized())

%timeit op()
%timeit nin17_for()
%timeit nin17_vectorized()


输出量:

5682372.071238058
5682372.071238058
5682372.071238059
2.54 s ± 8.79 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
5.62 ms ± 87.9 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
1.6 ms ± 11.2 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

相关问题