Mudanças entre as edições de "Números Aleatórios"

De Física Computacional
Ir para: navegação, pesquisa
(Criou página com '[http://davinci.if.ufrgs.br/wiki/index.php/Números_aleatórios] Números Aleatórios')
 
 
(18 revisões intermediárias por 2 usuários não estão sendo mostradas)
Linha 1: Linha 1:
[http://davinci.if.ufrgs.br/wiki/index.php/Números_aleatórios] 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:
 +
<pre>
 +
# 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)
 +
</pre>
 +
 
 +
== 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:
 +
 
 +
<pre>
 +
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}
 +
</pre>
 +
 
 +
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.
 +
 
 +
<pre>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))
 +
</pre>
 +
 
 +
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:
 +
 
 +
<pre>
 +
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)     
 +
</pre>
 +
 
 +
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.
 +
 
 +
<pre>
 +
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()
 +
</pre>
 +
 
 +
== Gás na caixa ==
 +
E por fim, podemos implementar a simulação do gás na caixa em Python da seguinte forma:
 +
<pre>
 +
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)
 +
</pre>
 +
 
 +
 
 +
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 <math display="inline">t</math> como <math display="inline">n_{+}\left(t\right)</math> então escrevemos:
 +
 
 +
<math display="block">n_{+}\left(t+1\right)=n_{+}\left(t\right)-p_{+}\left(t\right)+p_{-}\left(t\right)</math>
 +
 
 +
Onde <math display="inline">p_{+}\left(t\right)</math> é a probabilidade de encontrar uma partícula no lado direito no instante <math display="inline">t</math>. Ou seja, uma vez que qualquer partícula pode ser selecionada com a mesma probabilidade, então:
 +
 
 +
<math display="block">p_{+}\left(t\right)=\frac{n_{+}\left(t\right)}{N},\qquad p_{+}\left(t\right)=1-\frac{n_{+}\left(t\right)}{N},</math>
 +
 
 +
Substituindo:
 +
 
 +
<math display="block">\begin{align}
 +
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{align}</math>
 +
 
 +
Esta é a taxa de variação em um passo <math display="inline">\Delta t=1</math>. Então podemos obter o seguinte modelo de campo médio:
 +
 
 +
<math display="block">\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)</math>
 +
 
 +
<pre>
 +
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()               
 +
</pre>
 +
 
 +
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:
 +
 
 +
<math display="block">\left\langle n_{+}\right\rangle =\frac{1}{\tau}\sum_{t=t_{i}}^{t_{i}+\tau}n_{+}\left(t\right)</math>
 +
 
 +
Então a variância pode ser escrita como:
 +
 
 +
<math display="block">\sigma^{2}=\frac{1}{\tau}\sum_{t=t_{i}}^{t_{i}+\tau}\left(n_{+}\left(t\right)-\left\langle n_{+}\right\rangle \right)^{2}</math>
 +
 
 +
Uma vez que:
 +
 
 +
<math display="block">\begin{align}
 +
\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{align}</math>
 +
 
 +
Logo:
 +
 
 +
<math display="block">\sigma^{2}=\left\langle n_{+}^{^{2}}\left(t\right)\right\rangle -\left\langle n_{+}\right\rangle ^{2}</math>
 +
 
 +
== Distribuição binomial ==
 +
 
 +
Voltando ao exemplo do caminhante aleatório, considerando que tem <math display="inline">p</math> probabiliade de ir a direita (<math display="inline">q=1-p</math> de ir à esquerda), após <math display="inline">N</math> passos, tendo dado <math display="inline">n_{1}</math>passos à direita (<math display="inline">n_{2}=N-n_{1}</math> à esquerda), então a sua posição final é dada por <math display="inline">m=n_{1}-n_{2}</math>.
 +
 
 +
Como teve <math display="inline">N=n_{1}+n_{2}</math> passos, então sbustituindo em <math display="inline">m=n_{1}-n_{2}</math> podemos reescrever:
 +
 
 +
<math display="block">n_{1}=\frac{N+m}{2}\qquad n_{2}=\frac{N-m}{2}</math>
 +
 
 +
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 <math display="inline">n_{1}</math>pulos pra a direita e <math display="inline">n_{2}</math>pulos para esquerda é apenas:
 +
 
 +
<math display="block">p^{n_{1}}q^{n_{2}}=p^{\frac{N+m}{2}}q^{\frac{N-m}{2}}</math>
 +
 
 +
Agora precisamos multiplicar pelo número total de caminhos possíveis tendo <math display="inline">n_{1}</math>passos para a direita e <math display="inline">n_{2}</math>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 <math display="inline">m</math> após <math display="inline">N</math> passos, independente do caminho.
 +
 
 +
Isto é o número de modos que podemos colocar <math display="inline">n_{1}</math> objetos (de um total <math display="inline">N</math>) em uma ’caixa’ <math display="inline">n_{2}</math> em outra, lembrando que os objetos são essencialmente idênticos, só muda onde estamos guardando. Isto é de um conjunto de <math display="inline">N</math> “saltos” que vamos fazer, podemos escolher qualquer um dos <math display="inline">N</math> como a primeira escolha, então temos <math display="inline">N-1</math> opções como a segunda escolha, e sucessivamente, ao total temos <math display="inline">N!</math> arranjos possíveis. Porém todos os saltos à esquerda (ou bolinhas em uma caixa) são essencialmente os mesmos, e temos <math display="inline">n_{2}!</math> ordens possíveis de salto a esquerda, sendo que todas são iguais, do mesmo modo temos <math display="inline">n_{1}!</math> ordens possíveis de saltos a direita.
 +
 
 +
Para facilitar, vamos ilustrar com <math display="inline">N=4</math> e <math display="inline">n_{1}=n_{2}=2</math>. Isto é <math display="inline">m=0</math>. temos quatro saltos, 2 à esquerda <math display="inline">\left\{ s_{E}^{1},s_{E}^{2}\right\}</math> e 2 a direita <math display="inline">\left\{ S_{D}^{1},S_{D}^{2}\right\}</math>. Há <math display="inline">4!=24</math> formas de começar na posição <math display="inline">0</math> e terminar <math display="inline">0</math>, em um primeir momentos temos <math display="inline">4</math> saltos diferentes para escolher, depois vamos ter <math display="inline">3</math>, então <math display="inline">2</math> .... Por exemplo se escolhemos primeiro <math display="inline">s_{E}^{1}</math>, nos resta <math display="inline">\left\{ s_{E}^{2},S_{D}^{1},S_{D}^{2}\right\}</math>, que pode ser rearranjado de <math display="inline">6</math> formas diferentes:
 +
 
 +
<div class="center">
 +
 
 +
{|
 +
!align="center"| <math display="inline">\left\{ s_{E}^{1},s_{E}^{2},S_{D}^{1},S_{D}^{2}\right\}</math>
 +
|-
 +
|align="center"| <math display="inline">\left\{ s_{E}^{1},s_{E}^{2},S_{D}^{2},S_{D}^{1}\right\}</math>
 +
|-
 +
|align="center"| <math display="inline">\left\{ s_{E}^{1},S_{D}^{1}s_{E}^{2},S_{D}^{2}\right\}</math>
 +
|-
 +
|align="center"| <math display="inline">\left\{ s_{E}^{1},S_{D}^{1}s_{D}^{2},S_{E}^{2}\right\}</math>
 +
|-
 +
|align="center"| <math display="inline">\left\{ s_{E}^{1},S_{D}^{1},s_{E}^{2},S_{D}^{2}\right\}</math>
 +
|-
 +
|align="center"| <math display="inline">\left\{ s_{E}^{1},S_{D}^{1}S_{D}^{2},s_{E}^{2}\right\}</math>
 +
|}
 +
 
 +
 
 +
</div>
 +
O mesmo para cada um das outras <math display="inline">3</math> opções de inicício possível, então <math display="inline">4\times6=24</math>. Porém se a ordem de escolha dos saltos à direita for <math display="inline">\left\{ s_{D}^{1},s_{D}^{2}\right\}</math> ou <math display="inline">\left\{ s_{D}^{2},s_{D}^{1}\right\}</math>, não faz diferença. Então temos <math display="inline">2</math> ordens que são a indistinguíveis.
 +
 
 +
<div class="center">
 +
 
 +
{|
 +
!align="center"| <math display="inline">\left\{ s_{E}^{1},s_{E}^{2},S_{D},S_{D}\right\}</math>
 +
|-
 +
|align="center"| <math display="inline">\left\{ s_{E}^{1},S_{D}s_{E}^{2},S_{D}\right\}</math>
 +
|-
 +
|align="center"| <math display="inline">\left\{ s_{E}^{1},S_{D}s_{D},S_{E}^{2}\right\}</math>
 +
|}
 +
 
 +
 
 +
</div>
 +
Logo reduzimos pela metade as opções <math display="inline">\frac{4!}{2!}</math>. E o mesmo ainda pode ser dito para os saltos à esquerda. De forma que temos apenas <math display="inline">\frac{4!}{2!2!}=\frac{24}{4}=6</math> opções possíveis, que são:
 +
 
 +
<div class="center">
 +
 
 +
{|
 +
!align="center"| <math display="inline">\left\{ s_{E},s_{E},S_{D},S_{D}\right\}</math>
 +
|-
 +
|align="center"| <math display="inline">\left\{ s_{D},s_{D},S_{E},S_{E}\right\}</math>
 +
|-
 +
|align="center"| <math display="inline">\left\{ s_{E},s_{D},S_{E},S_{D}\right\}</math>
 +
|-
 +
|align="center"| <math display="inline">\left\{ s_{D},s_{E},S_{D},S_{E}\right\}</math>
 +
|-
 +
|align="center"| <math display="inline">\left\{ s_{E},s_{D},S_{D},S_{E}\right\}</math>
 +
|-
 +
|align="center"| <math display="inline">\left\{ s_{D},s_{E},S_{E},S_{D}\right\}</math>
 +
|}
 +
 
 +
 
 +
</div>
 +
Então o número total de caminhos indistinguiveis com <math display="inline">n_{1}</math>passos à direta e <math display="inline">n_{2}</math> à esquerda é: <math display="block">\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)!}</math>
 +
 
 +
Temos então que a probabilidade estar na posição <math display="inline">m</math> após <math display="inline">N</math> passos é dado por:
 +
 
 +
<math display="block">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}}</math>
 +
 
 +
Ou reescrevendo <math display="inline">n_{1}=n</math> para facilitar, e lembrando que <math display="inline">n_{1}=\frac{N+m}{2}</math> e <math display="inline">N=n_{1}+n_{2}</math>
 +
 
 +
<math display="block">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}</math>
 +
 
 +
Onde <math display="inline">\left(\begin{array}{c}
 +
N\\
 +
n
 +
\end{array}\right)=C_{n}^{N}</math> tabém é chamado de combinatória.
 +
 
 +
Utilizando o binômio de Newton:
 +
 
 +
<math display="block">\left(x+y\right)^{n}=\sum_{k=0}^{n}\left(\begin{array}{c}
 +
n\\
 +
k
 +
\end{array}\right)x^{n-k}y^{k}</math>
 +
 
 +
Temos a normalização:
 +
 
 +
<math display="block">\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</math>
 +
 
 +
Uma vez que <math display="inline">p+q=1</math>. E lembrando que <math display="inline">p\left(n\right)</math> é a probabilidade de estar em <math display="inline">m=2n-N</math>, ou seja de sair <math display="inline">n</math> vezes passo a direita em <math display="inline">N</math> tentativas, então o valor médio de <math display="inline">n</math> pode ser dado por:
 +
 
 +
<math display="block">\left\langle n\right\rangle =\sum_{n=0}^{N}np\left(n\right)</math>
 +
 
 +
E derivando o binômio de Newton:
 +
 
 +
<math display="block">\begin{array}{c}
 +
p\frac{d}{dp}\left(p+q\right)^{N}=pN\left(p+q\right)^{N-1}=pN\end{array}</math>
 +
 
 +
E também:
 +
 
 +
<math display="block">\begin{align}
 +
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{align}</math>
 +
 
 +
Então:
 +
 
 +
<math display="block">pN=\sum_{n=0}^{n}np\left(n\right)=\left\langle n\right\rangle</math>
 +
 
 +
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:
 +
 
 +
<math display="block">\sigma^{2}=\left\langle n^{2}\right\rangle -\left\langle n\right\rangle ^{2}</math>
 +
 
 +
Então:
 +
 
 +
<math display="block">\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}</math>
 +
 
 +
E:
 +
 
 +
<math display="block">\begin{align}
 +
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{align}</math>
 +
 
 +
Então:
 +
 
 +
<math display="block">p^{2}N^{2}-p^{2}N=\left\langle n^{2}\right\rangle -\left\langle n\right\rangle^2</math>
 +
 
 +
Ou ainda:
 +
 
 +
<math display="block">\begin{align}
 +
\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{align}</math>
 +
 
 +
Substituindo então:
 +
 
 +
<math display="block">\begin{align}
 +
\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{align}</math>
 +
 
 +
Logo:
 +
 
 +
<math display="block">\sigma=\sqrt{Npq}</math>
 +
 
 +
Além disso, para o limite <math display="inline">n\rightarrow\infty</math>, sendo <math> p \not\approx 1 </math>, temos a gaussiana <ref> [http://staff.ustc.edu.cn/~chenzyn/lectures/L.%20E.%20Reichl-A%20modern%20course%20in%20statistical%20physics.pdf A modern course in statistical physics] (L. E. Reichl)</ref>:
 +
 
 +
<math display="block">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)</math>
 +
 
 +
=== Exemplo: Caminhante aleatório ===
 +
 
 +
Voltando ao caminhante leatório, temos então que a posição média média no tempo <math display="inline">m\left(t\right)</math> vai ser <math display="inline">\left\langle m\right\rangle =0</math> se temos <math display="inline">p=q</math>, e se <math display="inline">p\neq q</math> então <math display="inline">\left\langle m\right\rangle =\left(p-q\right)N</math>. Além disso, a posição final é efetivamente <math display="inline">m=2n-N</math>, pois denotando<math display="inline">n_{D}</math> e <math display="inline">n_{E}</math> respectivamente como a quantidade de passos á direita e esquerda, temos:
 +
 
 +
<math display="block">n_{D}-n_{E}=n_{D}-\left(N-n_{D}\right)=2n_{D}-N</math>
 +
 
 +
Logo nossa média é:
 +
 
 +
<math display="block">\begin{align}
 +
\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{align}</math>
 +
 
 +
Logo:
 +
 
 +
<math display="block">\left\langle m\right\rangle ^{2}=4\left\langle n\right\rangle ^{2}-4\left\langle n\right\rangle N+N^{2}</math>
 +
 
 +
E:
 +
 
 +
<math display="block">\begin{align}
 +
\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{align}</math>
 +
 
 +
Fazeno então o desvio:
 +
 
 +
<math display="block">\begin{align}
 +
\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{align}</math>
 +
 
 +
Temos então que:
 +
 
 +
<math display="block">\sigma_{m}=2\sigma_{n}</math>
 +
 
 +
Ou simplesmente:
 +
 
 +
<math display="block">\sigma_{m}=2\sqrt{pqN}</math>
 +
 
 +
E sendo <math display="inline">p=q=\frac{1}{2}</math>:
 +
 
 +
<math display="block">\sigma_{m}=\sqrt{N}</math>
 +
 
 +
Sendo ainda <math display="inline">N</math> corresponde ao tempo, então fazendo <math display="inline">N=t</math>
 +
 
 +
<math display="block">\sigma_{m}=\sqrt{t}</math>
 +
 
 +
<pre>
 +
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()                   
 +
</pre>
 +
 
 +
== Distribuição de Poisson ==
 +
 
 +
Temos novamente <math display="inline">N\rightarrow\infty</math>, mas agora sendo <math display="inline">p\rightarrow0</math>, onde <math display="inline">pN=\lambda>0</math>. Tendo um evento que é continuamente testado durante um intervalo de tempo <math display="inline">\tau</math>, se a quantidade de testes neste intervalo é escrito como <math display="inline">N=\frac{\tau}{\delta t}</math>, tratando a probabilidade de maneira contínua <math display="inline">p\rightarrow\delta p</math> ficamos com:
 +
 
 +
<math display="block">\lambda=Np=\tau\frac{\delta p}{\delta t}</math>
 +
 
 +
<math display="inline">\frac{\delta p}{\delta t}</math> é probabilidade por tempo, e <math display="inline">\tau</math> o tempo total do intervalo, logo <math display="inline">\lambda</math> é o número médio de eventos no período <math display="inline">\tau</math>. Além disso, com <math display="inline">N=1</math>, temos o tempo médio de um evento <math display="inline">\tau=1/\lambda</math>. Retomando então a distribuição binomial:
 +
 
 +
<math display="block">p\left(n\right)=\left(\begin{array}{c}
 +
N\\
 +
n
 +
\end{array}\right)p^{n}q^{N-n}</math>
 +
 
 +
E manipulando:
 +
 
 +
<math display="block">\begin{align}
 +
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{align}</math>
 +
 
 +
No limite de <math display="inline">N\rightarrow\infty</math> então:
 +
 
 +
<math display="block">\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}</math>
 +
 
 +
E se <math display="inline">\lim_{N\rightarrow\infty}\left(N-a\right)\approx\lim_{N\rightarrow\infty}N</math> sendo <math display="inline">a</math> uma constante qualquer, então, como temos <math display="inline">n</math> termos no numerador e <math display="inline">N^{n}</math> nos deixa com <math display="inline">n</math> termos no denominador:
 +
 
 +
<math display="block">\begin{align}
 +
\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{align}</math>
 +
 
 +
Temos então:
 +
 
 +
<math display="block">p_{\lambda}\left(n\right)=\frac{e^{-\lambda}\lambda^{n}}{n!}</math>
 +
O valor mais provável é <math display="inline">\left\langle n\right\rangle =\lambda</math> e temos que <math display="inline">\sigma^{2}=\lambda</math>. Pois obtemos do resultado da binomial que <math display="inline">\left\langle n\right\rangle =pN</math> e <math display="inline">\sigma=\sqrt{pN\left(1-p\right)}</math>. No limite em que <math display="inline">p\rightarrow0</math> e <math display="inline">N\rightarrow \infty</math> temos então que <math display="inline">\left\langle n\right\rangle =pN=\lambda</math> e:
 +
 
 +
 
 +
<math display="block">\begin{align}
 +
\sigma= & \lim_{p\rightarrow0}\sqrt{pN\left(1-p\right)}=\lim_{p\rightarrow0}\sqrt{\lambda\left(1-p\right)}=\sqrt{\lambda}\end{align}</math>
 +
 
 +
Ou seja <math display="inline">\sigma^{2}=\lambda</math>
 +
 
 +
 
 +
 
 +
<pre>
 +
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()
 +
</pre>
 +
 
 +
= Citações =
 +
<references />

Edição atual tal como às 23h18min de 6 de maio de 2022

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 como então escrevemos:

Onde é a probabilidade de encontrar uma partícula no lado direito no instante . Ou seja, uma vez que qualquer partícula pode ser selecionada com a mesma probabilidade, então:

Substituindo:

Esta é a taxa de variação em um passo . Então podemos obter o seguinte modelo de campo médio:

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:

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

Uma vez que:

Logo:

Distribuição binomial

Voltando ao exemplo do caminhante aleatório, considerando que tem probabiliade de ir a direita ( de ir à esquerda), após passos, tendo dado passos à direita ( à esquerda), então a sua posição final é dada por .

Como teve passos, então sbustituindo em podemos reescrever:

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 pulos pra a direita e pulos para esquerda é apenas:

Agora precisamos multiplicar pelo número total de caminhos possíveis tendo passos para a direita e 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 após passos, independente do caminho.

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

Para facilitar, vamos ilustrar com e . Isto é . temos quatro saltos, 2 à esquerda e 2 a direita . Há formas de começar na posição e terminar , em um primeir momentos temos saltos diferentes para escolher, depois vamos ter , então .... Por exemplo se escolhemos primeiro , nos resta , que pode ser rearranjado de formas diferentes:


O mesmo para cada um das outras opções de inicício possível, então . Porém se a ordem de escolha dos saltos à direita for ou , não faz diferença. Então temos ordens que são a indistinguíveis.


Logo reduzimos pela metade as opções . E o mesmo ainda pode ser dito para os saltos à esquerda. De forma que temos apenas opções possíveis, que são:


Então o número total de caminhos indistinguiveis com passos à direta e à esquerda é:

Temos então que a probabilidade estar na posição após passos é dado por:

Ou reescrevendo para facilitar, e lembrando que e

Onde tabém é chamado de combinatória.

Utilizando o binômio de Newton:

Temos a normalização:

Uma vez que . E lembrando que é a probabilidade de estar em , ou seja de sair vezes passo a direita em tentativas, então o valor médio de pode ser dado por:

E derivando o binômio de Newton:

E também:

Então:

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:

Então:

E:

Então:

Ou ainda:

Substituindo então:

Logo:

Além disso, para o limite , sendo , temos a gaussiana [1]:

Exemplo: Caminhante aleatório

Voltando ao caminhante leatório, temos então que a posição média média no tempo vai ser se temos , e se então . Além disso, a posição final é efetivamente , pois denotando e respectivamente como a quantidade de passos á direita e esquerda, temos:

Logo nossa média é:

Logo:

E:

Fazeno então o desvio:

Temos então que:

Ou simplesmente:

E sendo :

Sendo ainda corresponde ao tempo, então fazendo

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 , mas agora sendo , onde . Tendo um evento que é continuamente testado durante um intervalo de tempo , se a quantidade de testes neste intervalo é escrito como , tratando a probabilidade de maneira contínua ficamos com:

é probabilidade por tempo, e o tempo total do intervalo, logo é o número médio de eventos no período . Além disso, com , temos o tempo médio de um evento . Retomando então a distribuição binomial:

E manipulando:

No limite de então:

E se sendo uma constante qualquer, então, como temos termos no numerador e nos deixa com termos no denominador:

Temos então:

O valor mais provável é e temos que . Pois obtemos do resultado da binomial que e . No limite em que e temos então que e:


Ou seja


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()

Citações

  1. A modern course in statistical physics (L. E. Reichl)