Fast Python/Numpy Frequency-Severity Distribution Simulation
<p>I'm looking for a away to simulate a classical frequency severity distribution, something like:
X = sum(i = 1..N, Y_i), where N is for example poisson distributed and Y lognormal.</p>

<p>Simple naive numpy script would be:</p>

<pre><code>import numpy as np
SIM_NUM = 3

X = []
for _i in range(SIM_NUM):
    nr_claims = np.random.poisson(1)
    temp = []
    for _j in range(nr_claims):
         temp.append(np.random.lognormal(0, 1))
    X.append(sum(temp))
</code></pre>

<p>Now I try to vectorize that for a performance increase. A bit better is the following:</p>

<pre><code>N = np.random.poisson(1, SIM_NUM)
X = []
for n in N:
    X.append(sum(np.random.lognormal(0, 1, n)))
</code></pre>

<p>My question is if I can still vectorize the second loop? For example by simulating all the losses with:</p>

<pre><code>N = np.random.poisson(1, SIM_NUM)
# print(N) would lead to something like [1 3 0]
losses = np.random.lognormal(0,1, sum(N))
# print(N) would lead to something like 
#[ 0.56750244  0.84161871  0.41567216  1.04311742]

# X should now be [ 0.56750244, 0.84161871 + 0.41567216 + 1.04311742, 0] 
</code></pre>

<p>Ideas that I have are "smart slicing" or "matrix multiplication with A = [[1, 0, 0, 0]],[0,1,1,1],[0,0,0,0]]. But I couldn't make something clever out these ideas. </p>

<p>I'm looking for fastest possible computation of X.</p>

<p>You're looking for <a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.ufunc.reduceat.html" rel="nofollow"><code>numpy.add.reduceat</code></a>:</p>

<pre><code>N = np.random.poisson(1, SIM_NUM)
losses = np.random.lognormal(0,1, np.sum(N))

x = np.zeros(SIM_NUM)
offsets = np.r_[0, np.cumsum(N[N&gt;0])]
x[N&gt;0] = np.add.reduceat(losses, offsets[:-1])
</code></pre>

<p>The case where <code>n == 0</code> is handled separately, because of how <code>reduceat</code> works. Also, be sure to use <code>numpy.sum</code> on arrays instead of the much slower Python <code>sum</code>.</p>

<p>If this is faster than the other answer depends on the mean of your Poisson distribution.</p>

python performance numpy vectorization