Números Aleatórios

Introdução

Existem na natureza processos de caráter estocástico, cujo resultado não pode ser predito a priori. Os exemplos típicos são o lançamento de uma moeda, o jogo dos dados ou a roleta. Em cada um desses casos sabemos quais são os resultados possíveis, porem não ha como predizer o resultado da próxima jogada. Entretanto podemos sim dizer qual será resultado em termos estatísticos de um grande número de jogadas: se lançarmos 1000 vezes uma moeda aproximadamente 50% das vezes o resultado será cara, e a percentagem será cada vez mais próxima de 50% quanto maior o número de lançamentos. O mesmo pode ser dito da roleta: num grande número de jogadas aproximadamente 1/37 das vezes terá dado o 17 por exemplo. Embora não seja suficiente para ganhar no jogo é um resultado previsível do ponto de vista estatístico. De uma outra forma, se nos fosse permitido colocar uma ficha nos pares e no zero ao mesmo tempo repetindo a jogada um grande número de vezes, teríamos lucro no final das contas.

Por outra parte existem na física processos cuja natureza é em principio determinística, cujo resultado é, em princípio, previsível conhecendo as condições inicias, mas que na prática resulta impossível como por exemplo a exata trajetória das moléculas de um gás.

Mesmo assim poderíamos integrar as equações de Newton de N partículas (moléculas ou planetas) com um dos algoritmos estáveis usados na unidade anterior: o Euler-Cromer ou de preferência o Verlet. O principio é o mesmo que no caso de uma única partícula como foi o caso do oscilador harmônico unidimensional. No caso de N partículas existem algumas complicações adicionais pois é necessário somar a força que as N-1 partículas fazem em uma dada partícula. À integração numérica das equações de movimento de N partículas é dada o nome de Dinâmica Molecular. O método é muito usado é viável num PC até 100000 ou 1 milhão de partículas. Porem depois de um determinado tempo a trajetória obtida para uma dada partícula pode estar muito afastada da trajetória teórica. No fundo acontece algo parecido com o problema da moeda. Ele é em princípio um problema determinístico mas na prática resulta randômico. Se diz que são problema sensíveis as condições iniciais: uma pequena variação na posição ou velocidade inicial da moeda se traduz numa trajetória muito diferente depois de um certo tempo.

Mas no caso de um gás, no fim não interessa pois nada importa a posição exata de uma dada molécula, o que importa são as propriedades estatísticas do sistema com um todo. Fazendo o análogo com o problema da moeda, é saber que proporção de cara ou coroa teremos depois de N lançamentos. No gás as medidas que interessam são a pressão, energia, temperatura, etc, propriedades médias do sistema em equilíbrio.

Dependendo das propriedades de interesse e dada essa caraterística randômica do movimento individual complexo de uma molécula qualquer, podemos adotar um modelo mais simples para simular o comportamento do sistema que o método de Dinâmica Molecular, relativamente difícil de programar. Tomemos o exemplo de um gás que está inicialmente contido numa metade de uma caixa, que está dividida ao meio por uma parede com um fim pequeno furo pelo qual pode passar uma molécula por vez. Sabemos que depois de um tempo o sistema alcançará o equilíbrio e que o número médio de partículas em cada metade da caixa será N/2.

Uma maneira simples de simular esse processo e supor que as partículas não interagem entre elas e que a probabilidade (por unidade de tempo) de uma dada partícula trocar de lado e a mesma para todas independente do lado da caixa na qual se encontre. Podemos implementar esse modelo escolhendo aleatoriamente uma partícula e trocá-la de lado. Acompanhando este processo no tempo deveríamos observar a evolução do sistema ao equilíbrio.

Números aleatórios

Para implementar esse modelo precisamos números aleatórios. Os verdadeiros números aleatórios são aqueles que resultam de um experimento estocástico, como os exemplos citados ao começo: moeda, dados, roleta. Mas para poder usar em cálculos no computador devemos achar um método mais rápido de gerar um seqüência de números aleatórios, deve ser o própio computador quem forneça essa seqüência. Isso é possível como veremos em seguida, porem eles não são estritamente aleatórios, eles parecem não obedecer a nenhum padrão, nenhuma seqüência lógica parece estar por trás desses números, mas na verdade se trata de uma seqüência completamente previsível, só para quem conhece a regra de geração da seqüência, ou seja ninguém poderia acertar o próximo número a menos que conheça a receita. Ou seja, que podem se comportar como números aleatórios tão bons quanto os verdadeiros. Por isso são as vezes chamados de pseudo-aleatórios.

Uma maneira de gerar uma seqüência assim é por meio da operação de congruência (resto da divisão inteira):

xn+1 = resto { (a xn + b)/m }

onde a, b, m e os xn são todos números inteiros. Dependendo da escolha dos inteiros a, b e m, a sequencia de números {x1, x2, ..., xn} é uma sequencia randómica dos números entre 1 e m. A melhor maneira de entender como isso funciona é com um exemplo: Seja a=4, b=3 e m=1, e seja o x inicial x0= 2

Em python, o método da congruência pode ser implementado da seguinte forma:

# gerador de num aleatorios de congruencia
def randi():
    global x
    a = 1029; b= 221591; m = 1048576
    x = (a*x + b) % m
    return x

# entrar aqui um interiro como semente inicial
x = int(input('semente (entrar inteiro)? '))
for i in range(10):
    n = randi()
    print(i,n)

Exemplos

Como uma solução nativa em Python temos a biblioteca random, para obter um número aleatório entre 0 e 1 podemos utilizar random(), ou randint(a,b) para inteiros entre (incluindo) a e b:

import random
print("Número aleatório: ",random.random()) # Frações no conjunto [0,1)
print("Roleta:",random.randint(0,37))       # Inteiros no conjunto [0,37]
print("Dado:",int(random.random()*6)+1)     # Inteiros no conjunto [1,6]
print("Moeda: ", 2*random.randint(0,1)-1)   # Inteiros no conjunto {-1,1}

O caminhante aleatório, pode ser facilmente implementado, tanto para se ter a mesma probabilidade de se mover em qualquer direção quanto para probabilidades diferentes.

import random
import matplotlib.pyplot as plt

n = 1000000 #Número de passos
p =0.5    #probabilidade de ir para o lado direito
x=[0]     #Posição inicial
c=0       #contador
for j in range(n):                #Com n passos em cada repetição
  dx = (+1) if (random.random()<=p) else (-1)
  x.append(x[j]+dx)               #Atualizamos a posição
  c = (c+1) if (dx>0) else (c)    #Atualizamos o contador
plt.plot(x)       
print("A particula se moveu para a direita {:.2f}% das vezes".format(100*c/n))

Também podemos implementar uma "apuração" de uma eleição. Por exemplo supondo que determinado candidato tem p% de chance de receber um voto, então:

import random
import matplotlib.pyplot as plt

p =0.5    #probabilidade do candidato ser votado
v=[0]     #Quantidade de votos em %
T=0       #Quantidade de votos totais
for i in range(1,1601): #Quantidade de candidatos
  x=(1) if (random.random()<=p) else (0)
  T+=x
  v.append(T/i)
plt.plot(v)       

Um histograma é um meio interessante para analisar alguns resultados. Considerando um jogo no qual uma moeda é jogada M vezes, e se repetimos o jogo por N vezes, podemo montar um histograma sobre a quantidade de caras que obtemos.

from matplotlib.ticker import FormatStrFormatter
import random
import matplotlib.pyplot as plt

N = 10000 #Quantidade de jogos
M = 100   #Quantidade de vezes que jogamos a moeda em cada jogo
res=[]    #Onde vamos guardar cada resultado
for i in range(N): #Vamos jogar N 
  s=0
  for j in range (M): #M moedas
    p=random.random()<0.5
    s = (s+1) if (p) else (s) #Soamos se saiu cara
  res.append(s) 
plt.hist(res,bins=[i for i in range(M+1)])
plt.show()

Gás na caixa

E por fim, podemos implementar a simulação do gás na caixa em Python da seguinte forma:

import random
import matplotlib.pyplot as plt

N = 1000     #Quantidade de partículas
cx =N*[0]    #Consideramos 0 à esquerda, e 1 à direita
np = 10000   #Número de pasos
t = [sum(cx)]#Quantidade de partícula a direita

for i in range(np):
  j=random.randint(0,N-1) #Pegamos um índice
  cx[j] = (1) if (cx[j]==0) else (0)
  t.append(sum(cx))

plt.plot(t)

A continução temos um exemplo de código FORTRAN também implementando a simulação do exemplo do gas na caixa. No exemplo simulamos 1000 partículas e fica como exercício decidir que variáveis imprimir e como calculá-las. Também se é necessário imprimir a cada passo da simulação.

!Aula5: Programa Gas, simulação com 1000 partículas
Program Gas
Implicit None
Real*8 Randi
Integer N, i, j, t, x
Parameter (N = 1000)
Integer Box(N)

Read(*,*) t, x  !t:tempo de simulação, x:semente de Randi

Do i = 1, N
   Box(i) = -1
End Do

Do j = 1, t
   i = N*Randi(x)+1
   Box(i) = -Box(i)
   Print*, j, 'deixo/vc(exercício 2)'
End Do

End Program Gas

Observações sobre FORTRAN:

  1. Parameter é para assignar (em tempo de compilação) um valor fixo a uma variável que não pode mudar durante a execução do programa, evitando o erro decorrente de uma alteração involuntária por distração, etc.
  2. O segundo novo elemento é o vector A(10) ou B(5,20).
     No caso que nos ocupa é Box(N); a utilidade dele fica evidente no própio programa: de uma vez só definimos uma variável a qual, por meio do indice, permite guardar o estado de um grande número de variáveis. Neste caso são as posições das 1000 partículas.

Exercícios:

  1. Graficar a evolução temporal das partículas a esquerda e a direita da parede (ver figura acima).
  2. Observe o "tempo" necessário para entrar em regime estacionário. Como varia esse tempo se mudarmos o número de partículas. Você pode estabelecer uma relação entre t e N?
  3. Implemente uma maneira de calcular a média de partículas num lado da caixa (precissa os dois?) e veja como varia no tempo.
  4. Depende do intervalo usado para calcular a média?
  5. Observe a variação no número de partículas no regime estacionário. Implemente uma forma de quantificar essa variação no seu programa. Ou seja a desvio cuadrático médio ou flutuação
  6. Quantifique a dependencia do desvio no regime estacionário com o número de partículas

Campo médio

Escrevendo então a quantidade de partículas no lado direto da caixa em um tempo ${\textstyle t}$ como ${\textstyle n_{+}\left(t\right)}$ então escrevemos:

$$n_{+}\left(t+1\right)=n_{+}\left(t\right)-p_{+}\left(t\right)+p_{-}\left(t\right)$$

Onde ${\textstyle p_{+}\left(t\right)}$ é a probabilidade de encontrar uma partícula no lado direito no instante ${\textstyle t}$. Ou seja, uma vez que qualquer partícula pode ser selecionada com a mesma probabilidade, então:

$$p_{+}\left(t\right)={\frac {n_{+}\left(t\right)}{N}},\qquad p_{+}\left(t\right)=1-{\frac {n_{+}\left(t\right)}{N}},$$

Substituindo:

$${\begin{aligned}n_{+}\left(t+1\right)&=n_{+}\left(t\right)-{\frac {n_{+}\left(t\right)}{N}}+1-{\frac {n_{+}\left(t\right)}{N}}\\&=n_{+}\left(t\right)+1-2{\frac {n_{+}\left(t\right)}{N}}\\n_{+}\left(t+1\right)-&n_{+}\left(t\right)=1-{\frac {2n_{+}\left(t\right)}{N}}\end{aligned}}$$

Esta é a taxa de variação em um passo ${\textstyle \Delta t=1}$. Então podemos obter o seguinte modelo de campo médio:

$${\frac {d}{dt}}n_{+}\left(t\right)=1-{\frac {2n_{+}\left(t\right)}{N}}\longrightarrow n_{+}\left(t\right)={\frac {N}{2}}\left(1-e^{-{\frac {2t}{N}}}\right)$$

import random
import matplotlib.pyplot as plt
import numpy as np

N = 10000     #Quantidade de partículas
cx =N*[0]     #Consideramos 0 à esquerda, e 1 à direita
al = [0]      #Solução analítica para a quantidade de partículas à direita
NP = 30000    #Número de pasos
t = [sum(cx)] #Quantidade de partícula a direita

for i in range(NP):
  j=random.randint(0,N-1)               #Pegamos um índice
  cx[j] = (1) if (cx[j]==0) else (0)    #Números aleatórios
  t.append(sum(cx))
  al.append((1-np.exp(-2*(i+1)/N))*N/2) #Campo médio

plt.plot(t,label="Número aleatórios")
plt.plot(al,label="Campo-médio")
plt.legend()                 

O campo médio nos dá a média de várias trajetórias no sistema, isto é, se somássemos todas trajetórias do sistema, o que obtemos representa a média de infinitas trajetórias. Logo o campo-médio não só representa a média de várias trajetórias com baixas partículas, mas a medida que aumentamos o número de partículas, ele se aproxima desta trajetória.

Valor médio e variância

Tendo o valor médio de partículas a direita dado por:

$$\left\langle n_{+}\right\rangle ={\frac {1}{\tau }}\sum _{t=t_{i}}^{t_{i}+\tau }n_{+}\left(t\right)$$

Então a variância pode ser escrita como:

$$\sigma ^{2}={\frac {1}{\tau }}\sum _{t=t_{i}}^{t_{i}+\tau }\left(n_{+}\left(t\right)-\left\langle n_{+}\right\rangle \right)^{2}$$

Uma vez que:

$${\begin{aligned}\sigma ^{2}&={\frac {1}{\tau }}\sum _{t=t_{i}}^{t_{i}+\tau }\left(n_{+}\left(t\right)-\left\langle n_{+}\right\rangle \right)^{2}\\&={\frac {1}{\tau }}\sum _{t=t_{i}}^{t_{i}+\tau }\left(n_{+}^{2}+\left\langle n_{+}\right\rangle ^{2}-2n_{+}\left\langle n_{+}\right\rangle \right)\\&={\frac {1}{\tau }}\sum _{t=t_{i}}^{t_{i}+\tau }n_{+}^{2}+{\frac {1}{\tau }}\sum _{t=t_{i}}^{t_{i}+\tau }\left\langle n_{+}\right\rangle ^{2}-{\frac {2}{\tau }}\sum _{t=t_{i}}^{t_{i}+\tau }n_{+}\left\langle n_{+}\right\rangle \\&=\left\langle n_{+}^{2}\right\rangle +{\frac {\left\langle n_{+}\right\rangle ^{2}}{\tau }}\sum _{t=t_{i}}^{t_{i}+\tau }\left(1\right)-2\left\langle n_{+}\right\rangle \left[{\frac {1}{\tau }}\sum _{t=t_{i}}^{t_{i}+\tau }n_{+}\right]\\&=\left\langle n_{+}^{2}\right\rangle +\left\langle n_{+}\right\rangle ^{2}{\frac {\tau }{\tau }}-2\left\langle n_{+}\right\rangle \left\langle n_{+}\right\rangle \end{aligned}}$$

Logo:

$$\sigma ^{2}=\left\langle n_{+}^{^{2}}\left(t\right)\right\rangle -\left\langle n_{+}\right\rangle ^{2}$$

Distribuição binomial

Voltando ao exemplo do caminhante aleatório, considerando que tem ${\textstyle p}$ probabiliade de ir a direita (${\textstyle q=1-p}$ de ir à esquerda), após ${\textstyle N}$ passos, tendo dado ${\textstyle n_{1}}$passos à direita (${\textstyle n_{2}=N-n_{1}}$ à esquerda), então a sua posição final é dada por ${\textstyle m=n_{1}-n_{2}}$.

Como teve ${\textstyle N=n_{1}+n_{2}}$ passos, então sbustituindo em ${\textstyle m=n_{1}-n_{2}}$ podemos reescrever:

$$n_{1}={\frac {N+m}{2}}\qquad n_{2}={\frac {N-m}{2}}$$

A probabilidade pra uma sequência de pulos é apenas o produto das probabilidades dos pulos indivíduais, pois são eventos independentes. Então a probabilidade de fazer ${\textstyle n_{1}}$pulos pra a direita e ${\textstyle n_{2}}$pulos para esquerda é apenas:

$$p^{n_{1}}q^{n_{2}}=p^{\frac {N+m}{2}}q^{\frac {N-m}{2}}$$

Agora precisamos multiplicar pelo número total de caminhos possíveis tendo ${\textstyle n_{1}}$passos para a direita e ${\textstyle n_{2}}$passos para a esquerda, pois na equação acima temos apenas a probabilidade de um caminho, e estamos interessados na probabilidade final de estar na posição ${\textstyle m}$ após ${\textstyle N}$ passos, independente do caminho.

Isto é o número de modos que podemos colocar ${\textstyle n_{1}}$ objetos (de um total ${\textstyle N}$) em uma ’caixa’ ${\textstyle n_{2}}$ em outra, lembrando que os objetos são essencialmente idênticos, só muda onde estamos guardando. Isto é de um conjunto de ${\textstyle N}$ “saltos” que vamos fazer, podemos escolher qualquer um dos ${\textstyle N}$ como a primeira escolha, então temos ${\textstyle N-1}$ opções como a segunda escolha, e sucessivamente, ao total temos ${\textstyle N!}$ arranjos possíveis. Porém todos os saltos à esquerda (ou bolinhas em uma caixa) são essencialmente os mesmos, e temos ${\textstyle n_{2}!}$ ordens possíveis de salto a esquerda, sendo que todas são iguais, do mesmo modo temos ${\textstyle n_{1}!}$ ordens possíveis de saltos a direita.

Para facilitar, vamos ilustrar com ${\textstyle N=4}$ e ${\textstyle n_{1}=n_{2}=2}$. Isto é ${\textstyle m=0}$. temos quatro saltos, 2 à esquerda ${\textstyle \left\{s_{E}^{1},s_{E}^{2}\right\}}$ e 2 a direita ${\textstyle \left\{S_{D}^{1},S_{D}^{2}\right\}}$. Há ${\textstyle 4!=24}$ formas de começar na posição ${\textstyle 0}$ e terminar ${\textstyle 0}$, em um primeir momentos temos ${\textstyle 4}$ saltos diferentes para escolher, depois vamos ter ${\textstyle 3}$, então ${\textstyle 2}$ .... Por exemplo se escolhemos primeiro ${\textstyle s_{E}^{1}}$, nos resta ${\textstyle \left\{s_{E}^{2},S_{D}^{1},S_{D}^{2}\right\}}$, que pode ser rearranjado de ${\textstyle 6}$ formas diferentes:

${\textstyle \left\{s_{E}^{1},s_{E}^{2},S_{D}^{1},S_{D}^{2}\right\}}$
${\textstyle \left\{s_{E}^{1},s_{E}^{2},S_{D}^{2},S_{D}^{1}\right\}}$
${\textstyle \left\{s_{E}^{1},S_{D}^{1}s_{E}^{2},S_{D}^{2}\right\}}$
${\textstyle \left\{s_{E}^{1},S_{D}^{1}s_{D}^{2},S_{E}^{2}\right\}}$
${\textstyle \left\{s_{E}^{1},S_{D}^{1},s_{E}^{2},S_{D}^{2}\right\}}$
${\textstyle \left\{s_{E}^{1},S_{D}^{1}S_{D}^{2},s_{E}^{2}\right\}}$

O mesmo para cada um das outras ${\textstyle 3}$ opções de inicício possível, então ${\textstyle 4\times 6=24}$. Porém se a ordem de escolha dos saltos à direita for ${\textstyle \left\{s_{D}^{1},s_{D}^{2}\right\}}$ ou ${\textstyle \left\{s_{D}^{2},s_{D}^{1}\right\}}$, não faz diferença. Então temos ${\textstyle 2}$ ordens que são a indistinguíveis.

${\textstyle \left\{s_{E}^{1},s_{E}^{2},S_{D},S_{D}\right\}}$
${\textstyle \left\{s_{E}^{1},S_{D}s_{E}^{2},S_{D}\right\}}$
${\textstyle \left\{s_{E}^{1},S_{D}s_{D},S_{E}^{2}\right\}}$

Logo reduzimos pela metade as opções ${\textstyle {\frac {4!}{2!}}}$. E o mesmo ainda pode ser dito para os saltos à esquerda. De forma que temos apenas ${\textstyle {\frac {4!}{2!2!}}={\frac {24}{4}}=6}$ opções possíveis, que são:

${\textstyle \left\{s_{E},s_{E},S_{D},S_{D}\right\}}$
${\textstyle \left\{s_{D},s_{D},S_{E},S_{E}\right\}}$
${\textstyle \left\{s_{E},s_{D},S_{E},S_{D}\right\}}$
${\textstyle \left\{s_{D},s_{E},S_{D},S_{E}\right\}}$
${\textstyle \left\{s_{E},s_{D},S_{D},S_{E}\right\}}$
${\textstyle \left\{s_{D},s_{E},S_{E},S_{D}\right\}}$

Então o número total de caminhos indistinguiveis com ${\textstyle n_{1}}$passos à direta e ${\textstyle n_{2}}$ à esquerda é:

$${\frac {N!}{n_{1}!n_{2}!}}={\frac {N!}{n_{1}!\left(N-n_{1}\right)!}}={\frac {N!}{\left({\frac {N+m}{2}}\right)!\left({\frac {N-m}{2}}\right)!}}$$

Temos então que a probabilidade estar na posição ${\textstyle m}$ após ${\textstyle N}$ passos é dado por:

$$p\left(m,N\right)={\frac {N!}{\left({\frac {N+m}{2}}\right)!\left({\frac {N-m}{2}}\right)!}}p^{\frac {N+m}{2}}q^{\frac {N-m}{2}}$$

Ou reescrevendo ${\textstyle n_{1}=n}$ para facilitar, e lembrando que ${\textstyle n_{1}={\frac {N+m}{2}}}$ e ${\textstyle N=n_{1}+n_{2}}$

$$p\left(n\right)={\frac {N!}{n!\left(N-n\right)!}}p^{n}q^{N-n}=\left({\begin{array}{c}N\\n\end{array}}\right)p^{n}q^{N-n}$$

Onde ${\textstyle \left({\begin{array}{c}N\\n\end{array}}\right)=C_{n}^{N}}$ tabém é chamado de combinatória.

Utilizando o binômio de Newton:

$$\left(x+y\right)^{n}=\sum _{k=0}^{n}\left({\begin{array}{c}n\\k\end{array}}\right)x^{n-k}y^{k}$$

Temos a normalização:

$$\sum _{n=0}^{n}p\left(n\right)=\sum _{n=0}^{n}\left({\begin{array}{c}N\\n\end{array}}\right)q^{N-n}p^{n}=\left(p+q\right)^{N}=1$$

Uma vez que ${\textstyle p+q=1}$. E lembrando que ${\textstyle p\left(n\right)}$ é a probabilidade de estar em ${\textstyle m=2n-N}$, ou seja de sair ${\textstyle n}$ vezes passo a direita em ${\textstyle N}$ tentativas, então o valor médio de ${\textstyle n}$ pode ser dado por:

$$\left\langle n\right\rangle =\sum _{n=0}^{N}np\left(n\right)$$

E derivando o binômio de Newton:

$${\begin{array}{c}p{\frac {d}{dp}}\left(p+q\right)^{N}=pN\left(p+q\right)^{N-1}=pN\end{array}}$$

E também:

$${\begin{aligned}p{\frac {d}{dp}}\left(p+q\right)^{N}&=p{\frac {d}{dp}}\left[\sum _{n=0}^{n}\left({\begin{array}{c}N\\n\end{array}}\right)q^{N-n}p^{n}\right]\\&=p\sum _{n=0}^{n}n\left({\begin{array}{c}N\\n\end{array}}\right)q^{N-n}p^{n-1}\\&=\sum _{n=0}^{n}n\left({\begin{array}{c}N\\n\end{array}}\right)q^{N-n}p^{n}\\&=\sum _{n=0}^{n}np\left(n\right)\end{aligned}}$$

Então:

$$pN=\sum _{n=0}^{n}np\left(n\right)=\left\langle n\right\rangle$$

Para o desvio padrão adotamos uma estratégia parecida mas com a segunda derivada do binômio de Newton. Sendo o desvio dado por:

$$\sigma ^{2}=\left\langle n^{2}\right\rangle -\left\langle n\right\rangle ^{2}$$

Então:

$${\begin{array}{c}p^{2}{\frac {d^{2}}{dp^{2}}}\left(p+q\right)^{N}=p^{2}N\left(N-1\right)\left(p+q\right)^{N-1}=p^{2}N\left(N-1\right)\end{array}}$$

E:

$${\begin{aligned}p^{2}{\frac {d^{2}}{dp^{2}}}\left(p+q\right)^{N}&=p^{2}{\frac {d^{2}}{dp^{2}}}\left[\sum _{n=0}^{n}\left({\begin{array}{c}N\\n\end{array}}\right)q^{N-n}p^{n}\right]\\&=p^{2}\sum _{n=0}^{n}n\left(n-1\right)\left({\begin{array}{c}N\\n\end{array}}\right)q^{N-n}p^{n-2}\\&=\sum _{n=0}^{n}\left(n^{2}-n\right)\left({\begin{array}{c}N\\n\end{array}}\right)q^{N-n}p^{n}\\&=\sum _{n=0}^{n}n^{2}\left({\begin{array}{c}N\\n\end{array}}\right)q^{N-n}p^{n}-\sum _{n=0}^{n}n\left({\begin{array}{c}N\\n\end{array}}\right)q^{N-n}p^{n}\\&=\left\langle n^{2}\right\rangle -\left\langle n\right\rangle \end{aligned}}$$

Então:

$$p^{2}N^{2}-p^{2}N=\left\langle n^{2}\right\rangle -\left\langle n\right\rangle ^{2}$$

Ou ainda:

$${\begin{aligned}\left\langle n^{2}\right\rangle &=p^{2}N^{2}-p^{2}N+\left\langle n\right\rangle \\&=p^{2}N^{2}-p^{2}N+pN\end{aligned}}$$

Substituindo então:

$${\begin{aligned}\sigma ^{2}&=\left\langle n^{2}\right\rangle -\left\langle n\right\rangle ^{2}\\&=p^{2}N^{2}-p^{2}N+pN-\left(pN\right)^{2}\\&=pN-p^{2}N\\&=pN\left(1-p\right)\\&=pNq\end{aligned}}$$

Logo:

$$\sigma ={\sqrt {Npq}}$$

Além disso, para o limite ${\textstyle n\rightarrow \infty }$, sendo $p\not \approx 1$, temos a gaussiana :

$$p\left(n\right)=\left({\begin{array}{c}N\\n\end{array}}\right)p^{n}q^{N-n}\approx {\frac {1}{\sqrt {2\pi Npq}}}\exp \left({\frac {-\left(n-Np\right)^{2}}{2Npq}}\right)$$

Exemplo: Caminhante aleatório

Voltando ao caminhante leatório, temos então que a posição média média no tempo ${\textstyle m\left(t\right)}$ vai ser ${\textstyle \left\langle m\right\rangle =0}$ se temos ${\textstyle p=q}$, e se ${\textstyle p\neq q}$ então ${\textstyle \left\langle m\right\rangle =\left(p-q\right)N}$. Além disso, a posição final é efetivamente ${\textstyle m=2n-N}$, pois denotando${\textstyle n_{D}}$ e ${\textstyle n_{E}}$ respectivamente como a quantidade de passos á direita e esquerda, temos:

$$n_{D}-n_{E}=n_{D}-\left(N-n_{D}\right)=2n_{D}-N$$

Logo nossa média é:

$${\begin{aligned}\left\langle m\right\rangle &=\sum _{n=0}^{N}mp\left(n\right)\\&=\sum _{n=0}^{N}\left(2n-N\right)p\left(n\right)\\&=2\sum _{n=0}^{N}np\left(n\right)-\sum _{n=0}^{N}Np\left(n\right)\\&=2\left\langle n\right\rangle -N\sum _{n=0}^{N}p\left(n\right)\\&=2\left\langle n\right\rangle -N\end{aligned}}$$

Logo:

$$\left\langle m\right\rangle ^{2}=4\left\langle n\right\rangle ^{2}-4\left\langle n\right\rangle N+N^{2}$$

E:

$${\begin{aligned}\left\langle m^{2}\right\rangle &=\sum _{n=0}^{N}m^{2}p\left(n\right)\\&=\sum _{n=0}^{N}\left(2n-N\right)^{2}p\left(n\right)\\&=4\sum _{n=0}^{N}n^{2}p\left(n\right)-4N\sum _{n=0}^{N}np\left(n\right)+N^{2}\sum _{n=0}^{N}p\left(n\right)\\&=4\left\langle n^{2}\right\rangle -4N\left\langle n\right\rangle +N^{2}\end{aligned}}$$

Fazeno então o desvio:

$${\begin{aligned}\sigma _{m}^{2}&=\left\langle m^{2}\right\rangle -\left\langle m\right\rangle ^{2}\\&=\left(4\left\langle n^{2}\right\rangle -4N\left\langle n\right\rangle +N^{2}\right)-\left(4\left\langle n\right\rangle ^{2}-4\left\langle n\right\rangle N+N^{2}\right)\\&=4\left(\left\langle n^{2}\right\rangle -\left\langle n\right\rangle ^{2}\right)\\&=4\sigma _{n}^{2}\end{aligned}}$$

Temos então que:

$$\sigma _{m}=2\sigma _{n}$$

Ou simplesmente:

$$\sigma _{m}=2{\sqrt {pqN}}$$

E sendo ${\textstyle p=q={\frac {1}{2}}}$:

$$\sigma _{m}={\sqrt {N}}$$

Sendo ainda ${\textstyle N}$ corresponde ao tempo, então fazendo ${\textstyle N=t}$

$$\sigma _{m}={\sqrt {t}}$$

import random
import matplotlib.pyplot as plt

M = 5000 #Quantidade de simulações
N = 100  #Número de passos
p = 0.5  #probabilidade de ir para o lado direito
x = (N+1)*[0]
x2 = x.copy()
for i in range(M):
  y=[x[0]]
  for j in range(N):                #Com n passos em cada repetição
    dy = (+1) if (random.random()<=p) else (-1)
    y.append(y[j]+dy)
  x=np.add(x,y)               #Somando os x
  x2=np.add(x2,np.square(y))  #Somando os x²
xmed=x/M   # <x >
x2med=x2/M # <x²>
sig=np.sqrt(x2med-np.square(xmed)) #(<x²>-<x>²)^(1/2)
plt.plot(xmed,label="<x>")
plt.plot(sig,label='Sigma')
plt.legend()                     

Distribuição de Poisson

Temos novamente ${\textstyle N\rightarrow \infty }$, mas agora sendo ${\textstyle p\rightarrow 0}$, onde ${\textstyle pN=\lambda >0}$. Tendo um evento que é continuamente testado durante um intervalo de tempo ${\textstyle \tau }$, se a quantidade de testes neste intervalo é escrito como ${\textstyle N={\frac {\tau }{\delta t}}}$, tratando a probabilidade de maneira contínua ${\textstyle p\rightarrow \delta p}$ ficamos com:

$$\lambda =Np=\tau {\frac {\delta p}{\delta t}}$$

${\textstyle {\frac {\delta p}{\delta t}}}$ é probabilidade por tempo, e ${\textstyle \tau }$ o tempo total do intervalo, logo ${\textstyle \lambda }$ é o número médio de eventos no período ${\textstyle \tau }$. Além disso, com ${\textstyle N=1}$, temos o tempo médio de um evento ${\textstyle \tau =1/\lambda }$. Retomando então a distribuição binomial:

$$p\left(n\right)=\left({\begin{array}{c}N\\n\end{array}}\right)p^{n}q^{N-n}$$

E manipulando:

$${\begin{aligned}p\left(n\right)&={\frac {N!}{\left(N-n\right)!n!}}p^{n}q^{N-n}\\&={\frac {N!}{\left(N-n\right)!n!}}{\frac {\left(Np\right)^{n}}{N^{n}}}\left(1-{\frac {Np}{N}}\right)^{N-n}\\&=\left[{\frac {N!}{\left(N-n\right)!}}{\frac {1}{N^{n}n!}}\right]\lambda ^{n}\left(1-{\frac {\lambda }{N}}\right)^{N}\left(1-{\frac {\lambda }{N}}\right)^{-n}\\&=\left[{\frac {N\left(N-1\right)\left(N-2\right)\dots \left(N-n+1\right)\left(N-n\right)!}{\left(N-n\right)!}}\right]{\frac {\lambda ^{n}}{N^{n}n!}}\left(1-{\frac {\lambda }{N}}\right)^{N}\left(1-{\frac {\lambda }{N}}\right)^{-n}\\&={\frac {N\left(N-1\right)\dots \left(N-n+1\right)}{N^{n}}}{\frac {\lambda ^{n}}{n!}}\left(1-{\frac {\lambda }{N}}\right)^{N}\left(1-{\frac {\lambda }{N}}\right)^{-n}\end{aligned}}$$

No limite de ${\textstyle N\rightarrow \infty }$ então:

$$\lim _{N\rightarrow \infty }\left(1-{\frac {\lambda }{N}}\right)^{-n}=1,\qquad \lim _{N\rightarrow \infty }\left(1-{\frac {\lambda }{N}}\right)^{N}=e^{-\lambda }$$

E se ${\textstyle \lim _{N\rightarrow \infty }\left(N-a\right)\approx \lim _{N\rightarrow \infty }N}$ sendo ${\textstyle a}$ uma constante qualquer, então, como temos ${\textstyle n}$ termos no numerador e ${\textstyle N^{n}}$ nos deixa com ${\textstyle n}$ termos no denominador:

$${\begin{aligned}\lim _{N\rightarrow \infty }p\left(n\right)&=\lim _{N\rightarrow \infty }{\frac {N\left(N-1\right)\dots \left(N-n+1\right)}{N^{n}}}{\frac {\lambda ^{n}}{n!}}\left(1-{\frac {\lambda }{N}}\right)^{N}\left(1-{\frac {\lambda }{N}}\right)^{-n}\\&={\frac {\lambda ^{n}}{n!}}e^{-\lambda }\lim _{N\rightarrow \infty }\left[{\frac {N\left(N-1\right)\dots \left(N-n+1\right)}{N\dots N}}\right]\end{aligned}}$$

Temos então:

$$p_{\lambda }\left(n\right)={\frac {e^{-\lambda }\lambda ^{n}}{n!}}$$

O valor mais provável é ${\textstyle \left\langle n\right\rangle =\lambda }$ e temos que ${\textstyle \sigma ^{2}=\lambda }$. Pois obtemos do resultado da binomial que ${\textstyle \left\langle n\right\rangle =pN}$ e ${\textstyle \sigma ={\sqrt {pN\left(1-p\right)}}}$. No limite em que ${\textstyle p\rightarrow 0}$ e ${\textstyle N\rightarrow \infty }$ temos então que ${\textstyle \left\langle n\right\rangle =pN=\lambda }$ e:

$${\begin{aligned}\sigma =&\lim _{p\rightarrow 0}{\sqrt {pN\left(1-p\right)}}=\lim _{p\rightarrow 0}{\sqrt {\lambda \left(1-p\right)}}={\sqrt {\lambda }}\end{aligned}}$$

Ou seja ${\textstyle \sigma ^{2}=\lambda }$

import random
import matplotlib.pyplot as plt

N = 10000 #Quantidade de repetições
M = 1000  #Quantidade de passos em cada repetição
res=[]    #Onde vamos guardar cada resultado
lam = 3   #Lambda
p = lam/M #Probabilidade

for i in range(N): #Vamos jogar N 
  s=0
  for j in range (M): #M moedas
    b=random.random()<p
    s = (s+1) if (b) else (s) #Soamos se saiu cara
  res.append(s) 
plt.hist(res,bins=[j for j in range (lam-5,lam+6,1)])
plt.show()