perf: Use Python builtin random.random to generate Poisson samples

[rht]

I'm not sure of the typical values for the input to poisson_rvs.
There are 2 possible followup improvements:

1. For large mu, you can use the method described in https://stackoverflow.com/questions/69596814/improved-inverse-transform-method-for-poisson-random-variable-generation-in-r
2. The poisson_rvs could be defined in Cython instead.

Benchmark code:

import time
import random
import math

from scipy.stats import poisson

mean = 7.5

def poisson_rvs(mu):
    p0 = math.exp(-mu)
    F = p0
    i = 0
    sample = random.random()
    while sample >= F:
        i += 1
        F += p0 * (mu ** i) / math.factorial(i)
    return i


scipy_list = []
tic = time.time()
for i in range(1000):
    scipy_list.append(poisson.rvs(mean))
elapsed = time.time() - tic
print("scipy.stats", elapsed, sum(scipy_list) / len(scipy_list))
custom_list = []
tic = time.time()
for i in range(1000):
    custom_list.append(poisson_rvs(mean))
elapsed_custom = time.time() - tic
print("Custom", elapsed_custom, sum(custom_list) / len(custom_list))
print("Speedup", round(elapsed / elapsed_custom, 4))

scipy.stats 0.1656346321105957 7.525
Custom 0.008195161819458008 7.561
Speedup 20.2113

[rht]

Any update @nunezco2 ? I have one more followup perf PR, but waiting for this one to be merged first.

[snunezcr]

Hi rht, The change is very effective. Typical values for mu are not large enough to become problematic. I will perform a few other changes in the next weeks, and will address some of the issues with students in about a month. Regards,

…

On Mon, Jul 25, 2022 at 11:28 PM rht ***@***.***> wrote: I'm not sure of the typical values for the input to poisson_rvs. There are 2 improvements: 1. For large mu, you can use the method described in https://stackoverflow.com/questions/69596814/improved-inverse-transform-method-for-poisson-random-variable-generation-in-r 2. The poisson_rvs could be defined in Cython instead. Benchmark code: import timeimport randomimport math from scipy.stats import poisson mean = 7.5 def poisson_rvs(mu): p0 = math.exp(-mu) F = p0 i = 0 sample = random.random() while sample >= F: i += 1 F += p0 * (mu ** i) / math.factorial(i) return i scipy_list = []tic = time.time()for i in range(1000): scipy_list.append(poisson.rvs(mean))elapsed = time.time() - ticprint("scipy.stats", elapsed, sum(scipy_list) / len(scipy_list))custom_list = []tic = time.time()for i in range(1000): custom_list.append(poisson_rvs(mean))elapsed_custom = time.time() - ticprint("Custom", elapsed_custom, sum(custom_list) / len(custom_list))print("Speedup", round(elapsed / elapsed_custom, 4)) scipy.stats 0.1656346321105957 7.525 Custom 0.008195161819458008 7.561 Speedup 20.2113 ------------------------------ You can view, comment on, or merge this pull request online at: #37 Commit Summary - 87dc07b <87dc07b> perf: Use Python builtin random.random to generate Poisson samples File Changes (1 file <https://github.com/ncsa/COVID19-mesa/pull/37/files>) - *M* covidmodel.py <https://github.com/ncsa/COVID19-mesa/pull/37/files#diff-41ccf0e41c96b24529858505967abbcd46251c939334827f784ca6e6e8564b69> (16) Patch Links: - https://github.com/ncsa/COVID19-mesa/pull/37.patch - https://github.com/ncsa/COVID19-mesa/pull/37.diff — Reply to this email directly, view it on GitHub <#37>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AACHP2X2WOIG37QKPZV76MDVV5SPHANCNFSM54UPEBZQ> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

-- Santiago Núñez-Corrales, Ph.D. Urbana IL, USA www.linkedin.com/in/snunez +1 (217) 621-3470

[rht]

The current poisson_rvs is general. It's just that for large mu, there is an even better algorithm.