Random seed, comparing the output from two simple models

1

Here are two simple models that return different draws, even though they have the same random seed. The only difference is the indexing. I put PyTensor in fast compile mode in case it was changing the graph under the hood, to isolate that. It doesn’t change either output.

I sort of naively expect to get the same answers here, as you would get without the indexing. What might be happening under the hood? Or at least, I would expect to see the same values come out in the second model in a different order.

import pymc as pm
import numpy as np
import pytensor

random_seed = 42
mu = np.array([1,2,3,4,5])

with pytensor.config.change_flags(mode="FAST_COMPILE"):
    with pm.Model():
        x = pm.Normal("x", mu=mu, sigma=1.0)
        prior_samples = pm.draw(x, draws=3, random_seed=random_seed, mode="FAST_COMPILE")

prior_samples

array([[ 1.41832997,  2.60557617,  3.02878786,  2.915754  ,  6.46422098],
       [ 1.29072736,  0.66924358,  2.96527654,  4.28041847,  5.10749307],
       [-0.92080086,  3.57864499,  4.00595719,  4.45121505,  4.40656633]])

And then:

with pytensor.config.change_flags(mode="FAST_COMPILE"):
    with pm.Model():
        x = pm.Normal("x", mu=mu[:3], sigma=1.0)
        prior_samples = pm.draw(x, draws=3, random_seed=random_seed, mode="FAST_COMPILE")

prior_samples

array([[ 1.41832997,  2.60557617,  3.02878786],
       [-0.084246  ,  3.46422098,  3.29072736],
       [-0.33075642,  1.96527654,  3.28041847]])

You can see the first three elements of the first row are the same.

1 Like
2

You can get the same behavior with numpy. The random generator advances more times in the first function, so it ends up at a different place by the time the second call starts. The generator is not reset between draws

import numpy as np

mu = np.arange(5)

rng = np.random.default_rng(1)
draws1 = np.stack([rng.normal(mu) for _ in range(3)])
print(draws1)
# [[0.34558419 1.82161814 2.33043708 1.69684277 4.90535587]
#  [0.44637457 0.46304676 2.5811181  3.3645724  4.2941325 ]
#  [0.02842224 1.54671299 1.26354591 2.83709005 3.51788069]]

rng = np.random.default_rng(1)
draws2 = np.stack([rng.normal(mu[:3]) for _ in range(3)])
print(draws2)
# [[ 0.34558419  1.82161814  2.33043708]
#  [-1.30315723  1.90535587  2.44637457]
#  [-0.53695324  1.5811181   2.3645724 ]]

You can get any of these numbers manually like this (other distributions may not have a fixed number of rng steps):

nth_draw = 4

rng = np.random.default_rng(1)
rng.bit_generator.advance(nth_draw)
print(rng.normal(mu[nth_draw % 5]))
# 4.905355866673117

rng = np.random.default_rng(1)
rng.bit_generator.advance(nth_draw)
print(rng.normal(mu[:3][nth_draw % 3]))
# 1.9053558666731178
2 Likes