Numpy Random Choice with Non-regular Array Size

Question

I'm making an array of sums of random choices from a negative binomial distribution (nbd), with each sum being of non-regular length. Right now I implement it as follows:

import numpy
from numpy.random import default_rng
rng = default_rng()

nbd = rng.negative_binomial(1, 0.5, int(1e6))
gmc = [12, 35, 4, 67, 2]
n_pp = np.empty(len(gmc))
for i in range(len(gmc)):
    n_pp[i] = np.sum(rng.choice(nbd, gmc[i]))

This works, but when I perform it over my actual data it's very slow (gmc is of dimension 1e6), and I would like to vary this for multiple values of n and p in the nbd (in this example they're set to 1 and 0.5, respectively).

I'd like to work out a pythonic way to do this which eliminates the loop, but I'm not sure it's possible. I want to keep default_rng for the better random generation than the older way of doing it (np.random.choice), if possible.

Is there a reason you're sampling from a large, previously-generated sample each iteration instead of making a new negative_binomial call each time? — AJ Biffl
– AJ Biffl, Commented Nov 19, 2021 at 1:00
Yeah, as the answer below shows, I should be approaching this differently. I thought making a large sample base would speed things up as it's a one time generation, but this (obviously) doesn't leverage broadcasting. — HEP N008
– HEP N008, Commented Nov 19, 2021 at 4:00

Warren Weckesser · Accepted Answer · 2021-11-19 02:19:08Z

1

The distribution of the sum of m samples from the negative binomial distribution with parameters (n, p) is the negative binomial distribution with parameters (m*n, p). So instead of summing random selections from a large, precomputed sample of negative_binomial(1, 0.5), you can generate your result directly with negative_binomial(gmc, 0.5):

In [68]: gmc = [12, 35, 4, 67, 2]

In [69]: npp = rng.negative_binomial(gmc, 0.5)

In [70]: npp
Out[70]: array([ 9, 34,  1, 72,  7])

(The negative_binomial method will broadcast its inputs, so we can pass gmc as an argument to generate all the samples with one call.)

More generally, if you want to vary the n that is used to generate nbd, you would multiply that n by the corresponding element in gmc and pass the product to rng.negative_binomial.

edited Nov 19, 2021 at 2:19

answered Nov 19, 2021 at 1:42

Warren Weckesser

116k20 gold badges207 silver badges224 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

HEP N008 Over a year ago

I did it my unnecessary loop way and this way and compared histograms, and this works a treat. I'll have to look at the math underlying it a bit more to be sure I get it, but this appears to do the trick. Thanks!

Collectives™ on Stack Overflow

Numpy Random Choice with Non-regular Array Size

1 Answer 1

1 Comment

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

1 Comment

Your Answer

Sign up or log in

Post as a guest

Related