Outro teste às "regularidades"

Em vez de cronometrarmos a passagem de carros numa estrada (e ainda estou para o fazer, finalmente completando o artigo) ou coisa que o valha, podemos tentar testar as conclusões do artigo das regularidades através de uma simulação. É o que pretendo fazer agora.

Como simular computacionalmente uma série de instantes da ocorrência de um fenómeno aleatório num dado intervalo de tempo $I$? A maneira mais óbvia parece-me ser a seguinte: divide-se o intervalo de tempo dado num grande número de pequenos subintervalos e gera-se aleatoriamente a ocorrência (ou não) do fenómeno em cada um. Não me agrada este método porque envolve uma discretização temporal (que, na verdade, está sempre presente quando se usam computadores digitais, mas adiante), que introduz um parâmetro de escala (a duração dos subintervalos) que pode gerar efeitos irrelevantes.
Em vez disso, vamos gerar essa série simulada simplesmente ordenando por ordem crescente um conjunto de números aleatórios uniformemente distribuídos no intervalo $I$.
Em python, podemos fazê-lo com as instruções

import numpy as np
import matplotlib.pyplot as plt

tmin = 0
tmax = 3600
N = 36000
times = tmin + (tmax - tmin) * np.random.random(N)
times.sort()

(O módulo matplotlib será usado mais à frente, para gerar os gráficos.) O array times tem agora armazenados sequencialmente os instantes em que se dão 36000 ocorrências do fenómeno, no intervalo compreendido entre os instantes $t=0$ e $t=3600\,$s. Note que nesta simulação se dão 36000 ocorrências num intervalo com duração 3600 s, ou seja, em média dão-se 10 ocorrências por segundo. Como o número expectável de ocorrências de um fenómeno de Poisson com parâmetro $\lambda$ num intervalo de tempo com duração $t$ é $\lambda t$, concluímos que, nesta simulação, se toma $\lambda=10$ s^-1.

Probabilidade de $k$ ocorrências num intervalo $\delta t$

De acordo com a eq. (1) do artigo das regularidades, a probabilidade de que um fenómeno aleatório com contante de ocorrência $\lambda$ ocorra $k$ vezes num intervalo com duração $\delta t$ é \begin{equation} P_k(\delta t)=\frac{1}{k!}(\lambda\delta t)^k e^{-\lambda\delta t}. \end{equation} Vejamos se a nossa simulação satisfaz esta lei de probabilidade. Podemos estimar as probabilidades de $k$ ocorrências do fenómeno num intervalo de tempo $\delta t$ pela frequência relativa com que acontecem intervalos de duração $\delta t$ com $k$ ocorrências. OK, tomemos $\delta t=1s$. Vejamos quantas ocorrências ocorrem em cada segundo:

freqs = np.zeros(3600)
for i in times:
    freqs(int(t)) += 1

Este excerto começa por iniciar a zero um contador por cada segundo no intervalo $I$ (esses contadores são os elementos do array freqs); em seguida, corre todos os instantes armazenados no array times e, para cada um, incrementa o contador correspondente.
Contar agora quantos destes contadores registam zero, uma, duas, etc ocorrências consegue-se com a instrução

counts, numbers = np.unique(freqs,return_counts=True)

O array counts armazena os números de ocorrências por segundo e numbers os números de segundos (subintervalos) em que se dão esses números de ocorrências. Isto é, numbers[k] é o número de segundos na nossa simulação em que ocorrem counts[k] ocorrências. Calculemos agora as frequências relativas dos intervalos com $k$ ocorrências e as probabilidades dadas pela eq. (1) (temos que importar a função factorial do módulo scipy.misc):

from scipy.misc import factorial
cfreqs = numbers / 3600
probs = np.exp(-10)*(10)**counts/factorial(counts)

O fator 10 nesta última expressão é $\lambda\delta t$.
Falta agora representar graficamente estas variáveis. Nas instruções seguintes, constrói-se um gráfico em que as probabilidades $P_k(\delta t=1\,\text{s})$ são representadas por barras cinzentas e as frequências por cruzes azuis.

plt.xlabel(r'$k$'); plt.ylabel(r'$P_k(\delta t=1s)$');
plt.bar(counts,probs,color='0.75');
plt.plot(counts,freqs,'x');
plt.show()

Fora as flutuações expectáveis dada a natureza aleatória deste cálculo, o acordo entre os valores esperados e os obtidos é evidente. Fica assim verificado "experimentalmente" um resultado importante do artigo das regularidades. Outro resultado que poderia ser testado deste modo é o da densidade de probabilidade para a ocorrência de um hiato com duração $\delta t$. Fica como exercício para o leitor.