這是我想做的模擬：

1. 對於 1000 萬條模擬路徑中的每一條，
2. 我依次有 n = 100 個 lambda 值（所有路徑的 lambda 向量都相同），
3. 使用每個 lambda 值，我需要從Poisson分佈中生成一個隨機數，均值等於 lambda；
4. 假設Poisson分佈中的隨機數是$$ Poi(1), Poi(2), Poi(3)… Poi(100) $$,
5. 然後我需要從指數分佈中生成 100 組隨機數。在每組中，Poi(i) 是我需要的隨機數的個數。
6. 假設一組隨機數是： $$ Exp(1,1), Exp(1,2), …, Exp(1, Poi(1)) Exp(2,1), Exp(2,2), …, Exp(2, Poi(2)) … Exp(100,1), …, Exp(100,Poi(100)) $$
7. 每次模擬我需要的結果是一個數字，它是每組指數隨機數之和的加權平均值。即結果 = w(1)*$$ Exp(1,1)+Exp(1,2)+…+Exp(1,Poi(1)) $$+ … + 在 (100)$$ Exp(100,1)+Exp(100,2)+…+Exp(100,Poi(100)) $$

我怎樣才能以 Python 的方式做到這一點？我編寫了以下程式碼，它設法在每個 for 循環中進行一次模擬，每次模擬需要 1 秒。有沒有更快的模擬方法？

import numpy as np 
from math import pi 
import time

start = time.clock()

kau = 6.21 
theta = 0.019 
sigma_v = 0.61 
t = 1 
nSim = 10000000 
result = []

n = np.arange(1,7) 
lamb = 16 * pi**2 * n**2 / sigma_v**2 / t / (kau**2 * t**2 + 4 * pi**2 * n**2) 
gamma = (kau**2 * t**2 + 4 * pi**2 * n**2) / 2 / sigma_v **2 / t**2

np.random.seed(1)

for sim in range(nSim):
   Nn = np.random.poisson(lam=lamb)
   y = [np.random.exponential(1,size=i) for i in Nn]
   z = [sum(y[i]) for i in range(len(y))]
   X1 = sum(z/gamma)
   result.append(X1)

print("time elapsed = ",time.clock()-start)

如果我有時間，我稍後會對此進行破解，但基本上在 python 中最好的辦法是使用所有 RAM，例如，如果你想生成 10 毫米 RV，請不要這樣做：

for i in range(10e6):
   rv = stats.norm.rvs(size=1)

而是這樣做：

rvs = stats.norm.rvs(size=10e6)

避免 python 循環更可取，因此在使用之前存在變數會更好（如果 RAM 大小允許）。如果您的方程式/機制允許，您可以對所有計算進行矢量化處理，以便它們同時執行並返回所需值的數組。

你有這些循環：

for sim in range(nSim)和y = [np.random.exponential(1,size=i) for i in Nn] z = [sum(y[i]) for i in range(len(y))]，如果您使用適當變數的索引切片比固有循環可能要快得多（例如，即使您過度生成隨機變數並儲存在記憶體中並從中進行子選擇，它也會比在循環中精確生成您需要的變數更快）。

編輯：對於步驟 2-7

這是一個如何創建遮罩以子選擇頂部的範例 $ n_i $ 隨機數組的元素並將它們相加：

>>> np.random.seed(1)
>>> rvs = np.random.rand(16).reshape(4,4)
array([[4.17022005e-01, 7.20324493e-01, 1.14374817e-04, 3.02332573e-01],
      [1.46755891e-01, 9.23385948e-02, 1.86260211e-01, 3.45560727e-01],
      [3.96767474e-01, 5.38816734e-01, 4.19194514e-01, 6.85219500e-01],
      [2.04452250e-01, 8.78117436e-01, 2.73875932e-02, 6.70467510e-01]])
>>> index = np.tile(np.arange(1, 5)[:, np.newaxis], (1,4))
array([[1., 1., 1., 1.],
      [2., 2., 2., 2.],
      [3., 3., 3., 3.],
      [4., 4., 4., 4.]])
>>> number_in_col = np.array([1,2,1,3])
>>> mask = number_in_col >= index
array([[ True,  True,  True,  True],
      [False,  True, False,  True],
      [False, False, False,  True],
      [False, False, False, False]])
>>> e_sum = np.sum(rvs * mask, axis=0)
array([4.17022005e-01, 8.12663088e-01, 1.14374817e-04, 1.33311280e+00])
>>> weights = np.array([1.1, 1.2, 1.3, 2.2])
>>> simulation_value = np.sum(e_sum * weights)
4.36691675844615

上面的程式碼本質上是步驟 2-7 的框架。

要執行第 1 步，即重複 10 毫米次，您應該考慮 RAM 使用情況，並在數組中添加第三個軸，其中第三個軸包含有關每個模擬的資訊。您將獲得一維模擬值數組，而不是返回單個模擬值8.112136084，具體取決於您可以在一次通過中召集多少。

例如，具有 100 列的index(8-byte) rvs(8-byte) 和mask(1-byte) 記憶體，我猜 50 行佔每個模擬的 (1710050) 85,000 字節左右。如果您最多有 10GB RAM 空閒，您可以一次進行 100,000 次模擬並寫入數據，然後循環 100 次：

output = np.empty(shape=(100,100000))    
for i in range(100):
   # edit the above for third axis with 100,000 simulations each run
   output[i, :] = simulation_value[:]
return output.reshape(-1,)  # <- output is a 1D array of 10mm values.

我不打算測試它，但我保證它會比你目前的實現快得多。

編輯以擴展第 3 軸以獲得第 1 步的幫助

>>> rvs_ext = np.tile(rvs[:,:,np.newaxis], (1,1,3)) # <- create 3 copies along 3rd axis
>>> mask_ext = np.tile(mask[:,:,np.newaxis], (1,1,3))
>>> e_sum_ext = np.sum(rvs_ext * mask_ext, axis=0)
array([[4.17022005e-01, 4.17022005e-01, 4.17022005e-01],
      [8.12663088e-01, 8.12663088e-01, 8.12663088e-01],
      [1.14374817e-04, 1.14374817e-04, 1.14374817e-04],
      [1.33311280e+00, 1.33311280e+00, 1.33311280e+00]])
>>> simulation_value_ext = np.einsum('ij,i->j', e_sum_ext, weights)
array([4.36691676, 4.36691676, 4.36691676])

自從您平鋪原始值以來，您在這裡獲得了 3 個重複的模擬值，rvs但如果您嘗試了：

>>> rvs_ext = np.random.rand(48).reshape(4,4,3)

然後你會有不同的 sim 值。

編輯以考慮隨機Poisson向量。

要調整隨機 lambda 值，請嘗試以下操作：

>>> lambda_ext = np.array([[1,2,1,3], [2,2,3,3], [1,1,4,4]])
>>> lambda_3D = np.tile(lambda_ext.T[np.newaxis,:,:], (4,1,1)) 
>>> index_ext = np.tile(index[:,:,np.newaxis], (1,1,3))
>>> mask_3D = index_ext <= lambda_match
>>> mask_3D[:,:,0]
array([[ True,  True,  True,  True],
      [False,  True, False,  True],
      [False, False, False,  True],
      [False, False, False, False]])
>>> mask_3D[:,:,2]
array([[ True,  True,  True,  True],
      [False, False,  True,  True],
      [False, False,  True,  True],
      [False, False,  True,  True]])
>>> e_sum_ext = np.sum(rvs_ext * mask_3D, axis=0)

現在你有所有的骨架來組合一個腳本來執行這個。

引用自：https://quant.stackexchange.com/questions/41086