Aleatoriedade avançada

Fluxos independentes com `SeedSequence`

Quando várias partes do programa (ou vários processos) precisam de números aleatórios, não use sementes "consecutivas" (1, 2, 3...). O jeito documentado é partir de uma semente raiz e derivar fluxos independentes com spawn. A mesma raiz reproduz os mesmos fluxos, e eles não se sobrepõem:

import numpy as np

raiz = np.random.SeedSequence(2026)
geradores = [np.random.default_rng(filho) for filho in raiz.spawn(3)]
amostras = [g.random(2000) for g in geradores]
print(len(amostras), abs(np.corrcoef(amostras)[0, 1]) < 0.1)

a = np.random.default_rng(np.random.SeedSequence(2026).spawn(3)[1]).random(2)
b = np.random.default_rng(np.random.SeedSequence(2026).spawn(3)[1]).random(2)
print(np.array_equal(a, b))

3 True
True

O próprio gerador também sabe se dividir, com o método spawn, que devolve novos geradores independentes:

filhos = np.random.default_rng(1).spawn(2)
print(type(filhos[0]).__name__)

Generator

Em um programa paralelo, você cria os filhos no processo principal e entrega um a cada trabalhador. Assim cada um tem o seu fluxo, e o resultado total é reproduzível, não importa a ordem em que os trabalhadores terminem.

Guardar e restaurar o estado

O estado do gerador pode ser salvo e restaurado, o que permite retomar uma simulação do ponto em que parou:

rng = np.random.default_rng(1)
estado = rng.bit_generator.state
primeiro = rng.random(2)
rng.bit_generator.state = estado
print(np.array_equal(primeiro, rng.random(2)))

True

Simulação de Monte Carlo

Muitos problemas sem fórmula fechada se resolvem sorteando muitas vezes e contando. O clássico é estimar π: sorteie pontos em um quadrado de lado 1, e a fração que cai dentro do quarto de círculo vale π/4. A conta inteira é vetorizada:

import math


def estimar_pi(n, rng):
    xy = rng.random((n, 2))
    dentro = (xy ** 2).sum(axis=1) <= 1
    return 4 * dentro.mean()


estimativa = estimar_pi(1_000_000, np.random.default_rng(0))
print(abs(estimativa - math.pi) < 0.01)

True

O erro cai com a raiz quadrada do número de sorteios: para ganhar uma casa decimal de precisão, você precisa de cem vezes mais amostras. Dá para ver isso repetindo o experimento e medindo o quanto as estimativas variam:

def dispersao(n, repeticoes=30):
    rng = np.random.default_rng(3)
    return np.std([estimar_pi(n, rng) for _ in range(repeticoes)])


print("mais amostras, menos variação:", dispersao(100_000) < dispersao(1_000))

mais amostras, menos variação: True

Bootstrap: incerteza sem fórmula

Como saber o quanto confiar em uma média calculada de uma amostra? O bootstrap reamostra a própria amostra, com reposição, muitas vezes, calcula a média de cada reamostra e usa a dispersão delas como medida da incerteza. A matriz de reamostras é gerada de uma vez, e a média é tirada ao longo de um eixo:

rng = np.random.default_rng(5)
amostra = rng.normal(100, 15, 200)
medias = rng.choice(amostra, size=(2000, amostra.size), replace=True).mean(axis=1)
baixo, alto = np.percentile(medias, [2.5, 97.5])
print(baixo < amostra.mean() < alto, alto - baixo < 10)

True True

O intervalo entre os percentis 2,5% e 97,5% é um intervalo de confiança de 95% para a média. Aqui a largura ficou abaixo de 10, perto do esperado pela teoria (cerca de 4, pois o erro padrão é 15 / √200).

Amostrar de uma distribuição qualquer

Se você sabe a função de distribuição acumulada, a transformada inversa transforma sorteios uniformes em sorteios da distribuição que quiser. Para a exponencial, -ln(1 − u) / λ:

lam = 2.0
x = -np.log(1 - rng.random(100_000)) / lam
print(abs(x.mean() - 1 / lam) < 0.01)

True

Registre o gerador e a versão

Para um experimento que precisa ser repetido anos depois, guarde a semente raiz, o tipo de gerador (o padrão é o PCG64) e a versão do NumPy. A mesma semente reproduz os mesmos números na mesma versão, mas o projeto reserva o direito de mudar os algoritmos entre versões do Generator.