Question: I've got a problem at work that goes as follows: We have, say, 1 million possible events which define frequency-severity distributions. For each event we have an annual rate which defines a Poisson distribution, and alpha and beta parameters for a Beta distributions. The goal is to simulate in the order of >100,000 "years", each year being defined as getting a frequency N for each event, and getting N samples of the relative beta distribution. The cumbersome fact for me is how can I efficiently get N_i ~ Poisson(lambda_i) samples from Beta distribution Beta_i while also making sure I can attribute them to the correct year? In terms of output I'll need to look at both the maximum and total value of the samples per year, so temporarily I'm just storing it as an array of dictionaries (not intended to be the output format) years = 5000 rng = np.random.default_rng() losses = [] for year in range(years): occurences = rng.poisson(data['RATE']) annual_losses = [] for idx, occs in enumerate(occurences): if occs > 0: event = data.iloc[idx] for occ in range(occs): loss = rng.beta(event['alpha'], event['beta']) * event['ExpValue'] annual_losses.append(loss) annual_losses.append(0) losses.append({'year': year, 'losses': annual_losses})  I've tried to follow Optimising Python/Numpy code used for Simulation however I can't seem to understand how to vectorise this code efficiently. Changes I've made before posting here are (times for 5000 years):  swapping from scipy to numpy (72s -> 66s) calculating frequencies for all years in once outside the loop (66s -> 73s... oops)  Ideally I'd like this to run as fast, or for as many iterations, as possible, and I was also running into issues with memory using scipy previously. EDIT: Upon request here's a version that just calculates the maximum loss for each year years = 20_000 rng = np.random.default_rng() largest_losses = np.zeros(shape=years) for year in range(years): occurences = rng.poisson(data['RATE']) largest_loss = 0 for idx, event_occs in enumerate(occurences): if event_occs > 0: event = data.iloc[idx] for event_occ in range(event_occs): loss = rng.beta(event['alpha'], event['beta']) * event['ExpValue'] if loss > largest_loss: largest_loss = loss largest_losses[year] = largest_loss  For testing purposes the events have a total rate of ~0.997 and the timed tests above are from an event pool of 100,604 events. To specify the goal, I'd like to see if I can calculate the percentiles of the losses, e.g. the 1-in-250 loss np.percentile(largest_losses, (1 - 1/250) * 100)  which currently aren't accurate at 5,000, hence wanting a process that can run in a few minutes or less for ~100,000 years. EDIT 2: Reproduceable values (as np arrays): rate_scale = 4.81e-4 rate_shape = 2.06e-2 rates = rng.gamma(rate_shape, rate_scale, size=nevents) alpha_scale = 1.25 alpha_shape = 0.439 alphas = rng.gamma(alpha_shape, alpha_scale, size=nevents) beta_scale = 289 beta_shape = 0.346 betas = rng.gamma(beta_shape, beta_scale, size=nevents) value_scale = 4.2e8 value_shape = 1.0 values = rng.gamma(value_shape, value_scale, size=nevents)  and have amended my current code to use numpy arrays years = 20_000 rng = np.random.default_rng() largest_losses = np.zeros(shape=years) for year in range(years): occurences = rng.poisson(rates) largest_loss = 0 for idx, event_occs in enumerate(occurences): for event_occ in range(event_occs): loss = rng.beta(alphas[idx], betas[idx]) * values[idx] if loss > largest_loss: largest_loss = loss largest_losses[year] = largest_loss f'{np.percentile(largest_losses, (1 - 1/250) * 100):,.0f}' 

Answer: I profiled @WarrenWeckesser's answer, and I found that it spends 93% of its time generating poisson-distributed numbers using the rate array. Here's a line profile of what I'm talking about: Timer unit: 1e-09 s Total time: 3.57874 s File: /tmp/ipykernel_2774/1125349843.py Function: simulate_warren_answer at line 1 Line # Hits Time Per Hit % Time Line Contents ============================================================== 1 def simulate_warren_answer(nyears, rate, alpha, beta, value, rng): 2 1 13092.0 13092.0 0.0 max_loss = np.zeros(nyears, dtype=np.float64) 3 1 3150.0 3150.0 0.0 total_loss = np.zeros(nyears, dtype=np.float64) 4 1001 347152.0 346.8 0.0 for year in range(nyears): 5 1000 3355879564.0 3e+06 93.8 n = rng.poisson(rate) 6 1000 179118381.0 179118.4 5.0 k = np.nonzero(n)[0] 7 1000 911231.0 911.2 0.0 if len(k) > 0: 8 654 935876.0 1431.0 0.0 nonzero_n = n[k] 9 654 4049698.0 6192.2 0.1 a = np.repeat(alpha[k], nonzero_n) 10 654 2645836.0 4045.6 0.1 b = np.repeat(beta[k], nonzero_n) 11 654 2550944.0 3900.5 0.1 v = np.repeat(value[k], nonzero_n) 12 654 26183513.0 40036.0 0.7 losses = rng.beta(a, b)*v 13 654 3557136.0 5439.0 0.1 max_loss[year] = losses.max() 14 654 2543208.0 3888.7 0.1 total_loss[year] = losses.sum() 15 1 303.0 303.0 0.0 return max_loss, total_loss  For this reason, I started looking into whether it was possible to speed up poisson-distributed random number generation. I started by looking at the lengths of n and k, and realized that for the rates that the asker is interested in, k is much much smaller than n, and most entries in n are zero. If there is a 98% chance that the poisson-distributed number will be zero, I can simulate this by generating a uniformly distributed random number between [0, 1), and comparing it to 0.02. This leaves me to handle the other case 2% of the time, where the number will be 1+. I then realized that I could generalize this to generating any value of the poisson function, by comparing the random uniform value to the CDF. (In the code below, I am using the survival function, which is 1-CDF. This is for accuracy reasons; numbers near 0 are represented more accurately than numbers near 1.) This code generates a lookup table for converting uniformly-distributed numbers into poisson-distributed numbers according to each rate. Because the rates are small, most of the time, the code will only need to examine the first entry for each rate, and the result will be zero. It also does the k = np.nonzero(n)[0] step inside this loop, as it is more cache efficient. Because the rates are not changing, the lookup table only needs to be build once, and can be re-used between years. import scipy import numpy as np import numba as nb def simulate_new_rng(nyears, rate, alpha, beta, value, rng): max_loss = np.empty(nyears, dtype=np.float64) total_loss = np.zeros(nyears, dtype=np.float64) lookup = build_lookup(rate) for year in range(nyears): uniform = rng.random(len(rate)) n, nonzeros = uniform_to_poisson_plus_nonzero(uniform, lookup) n_nonzeros = n[nonzeros] a = np.repeat(alpha[nonzeros], n_nonzeros) b = np.repeat(beta[nonzeros], n_nonzeros) v = np.repeat(value[nonzeros], n_nonzeros) losses = rng.beta(a, b)*v max_loss[year] = losses.max(initial=0) total_loss[year] = losses.sum() return max_loss, total_loss def build_lookup(rate, ptol=1e-15): """Build uniform -> poisson lookup table based on survival function Does not include an entry for k if the probability of the most likely rate to be k is less than ptol.""" ret = [] for k in itertools.count(): sf = scipy.stats.poisson.sf(k, rate) # include first row of lookup table no matter what if sf.max() < ptol and k != 0: break ret.append(sf) return np.array(ret) @nb.njit() def uniform_to_poisson_plus_nonzero(numbers, lookup): assert numbers.shape[0] == lookup.shape[1] nonzeros = [] ret = np.empty(lookup.shape[1], dtype=np.int64) for j in range(lookup.shape[1]): number = numbers[j] for i in range(lookup.shape[0]): if number < lookup[i, j]: continue else: break else: raise Exception("failed to find match for survival function - ptol is too big") ret[j] = i if i != 0: nonzeros.append(j) return ret, np.array(nonzeros, dtype=np.int64)  Benchmark result: Original code: 24.4 s ± 308 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) Warrens answer: 3.44 s ± 18.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) My answer: 648 ms ± 17.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)  Appendix: Full benchmarking code import matplotlib.pyplot as plt import itertools import scipy import numba as nb import numpy as np import pandas as pd rng = np.random.default_rng() nevents = 100604 nyears = 1000 rate_scale = 4.81e-4 rate_shape = 2.06e-2 rates = rng.gamma(rate_shape, rate_scale, size=nevents) alpha_scale = 1.25 alpha_shape = 0.439 alphas = rng.gamma(alpha_shape, alpha_scale, size=nevents) beta_scale = 289 beta_shape = 0.346 betas = rng.gamma(beta_shape, beta_scale, size=nevents) value_scale = 4.2e8 value_shape = 1.0 values = rng.gamma(value_shape, value_scale, size=nevents) def simulate_orig_answer(nyears, rate, alpha, beta, value, rng): data = pd.DataFrame({'RATE': rate, 'alpha': alpha, 'beta': beta, 'ExpValue': value}) rng = np.random.default_rng() largest_losses = np.zeros(shape=nyears) total_loss = np.zeros(nyears, dtype=np.float64) for year in range(nyears): occurences = rng.poisson(data['RATE']) largest_loss = 0 sum_loss = 0 for idx, event_occs in enumerate(occurences): if event_occs > 0: event = data.iloc[idx] for event_occ in range(event_occs): loss = rng.beta(event['alpha'], event['beta']) * event['ExpValue'] if loss > largest_loss: largest_loss = loss sum_loss += loss largest_losses[year] = largest_loss total_loss[year] = sum_loss return largest_losses, total_loss def simulate_warren_answer(nyears, rate, alpha, beta, value, rng): max_loss = np.zeros(nyears, dtype=np.float64) total_loss = np.zeros(nyears, dtype=np.float64) for year in range(nyears): n = rng.poisson(rate) k = np.nonzero(n)[0] if len(k) > 0: nonzero_n = n[k] a = np.repeat(alpha[k], nonzero_n) b = np.repeat(beta[k], nonzero_n) v = np.repeat(value[k], nonzero_n) losses = rng.beta(a, b)*v max_loss[year] = losses.max() total_loss[year] = losses.sum() return max_loss, total_loss def simulate_new_rng(nyears, rate, alpha, beta, value, rng): max_loss = np.empty(nyears, dtype=np.float64) total_loss = np.zeros(nyears, dtype=np.float64) lookup = build_lookup(rate) for year in range(nyears): uniform = rng.random(len(rate)) n, nonzeros = uniform_to_poisson_plus_nonzero(uniform, lookup) n_nonzeros = n[nonzeros] a = np.repeat(alpha[nonzeros], n_nonzeros) b = np.repeat(beta[nonzeros], n_nonzeros) v = np.repeat(value[nonzeros], n_nonzeros) losses = rng.beta(a, b)*v max_loss[year] = losses.max(initial=0) total_loss[year] = losses.sum() return max_loss, total_loss def build_lookup(rate, ptol=1e-15): """Build uniform -> poisson lookup table based on survival function Does not include an entry for k if the probability of the most likely rate to be k is less than ptol.""" ret = [] for k in itertools.count(): sf = scipy.stats.poisson.sf(k, rate) # include first row of lookup table no matter what if sf.max() < ptol and k != 0: break ret.append(sf) return np.array(ret) @nb.njit() def uniform_to_poisson_plus_nonzero(numbers, lookup): assert numbers.shape[0] == lookup.shape[1] nonzeros = [] ret = np.empty(lookup.shape[1], dtype=np.int64) for j in range(lookup.shape[1]): number = numbers[j] for i in range(lookup.shape[0]): if number < lookup[i, j]: continue else: break else: raise Exception("failed to find match for survival function - ptol is too big") ret[j] = i if i != 0: nonzeros.append(j) return ret, np.array(nonzeros, dtype=np.int64) %timeit simulate_orig_answer(nyears, rates, alphas, betas, values, rng) %timeit simulate_warren_answer(nyears, rates, alphas, betas, values, rng) %timeit simulate_new_rng(nyears, rates, alphas, betas, values, rng)