aula de hoje:

    1. o método de Monte Carlo
    1. integração por Monte Carlo
    1. universo homogêneo e isotrópico: distribuição uniforme de pontos em uma esfera
    1. MCMC em uma dimensão

simulações adicionais para os interessados:

  1. distribuição exponencial e distribuição gama
  1. área por integração de Monte Carlo
  1. simulação por MCMC de uma distribuição exponencial
  1. random walk em duas dimensões

1. o método de Monte Carlo

  • objetivo: obter amostras de uma distribuição de probabilidades unidimensional \(P(x)\)
  • dado um número aleatório \(\gamma\) uniformemente distribuído entre 0 e 1, obter \(x\) \[ \gamma = \int_{-\infty}^x P(x^\prime) dx^\prime = F(x)\]

\(F(x)\): função cumulativa de \(P(x)\)

\(x(\gamma)\) é obtido resolvendo-se a equação \(F(x)-\gamma = 0\).

ex.: simulação de uma distribuição exponencial

para a distribuição exponencial (que é definida para \(x>0\)): \[ \gamma = \int_0^x e^{-x^\prime} dx^\prime = 1-e^{-x}\]
logo \[x = -\ln (1-\gamma)=-\ln \gamma^\prime\] onde \(\gamma^\prime\) é também um número aleatório uniformemente distribuído entre 0 e 1.

# para saber o tempo de processamento do script:
t00 = Sys.time()

# como faremos simulações, vamos adotar uma semente para a geração de números aleatórios para assegurar a reprodutibilidade dos resultados
set.seed(1234)
# vamos gerar 1000 pontos como acima (por default log é o logarítmo natural):
x=-log(runif(1000))

# estatísticas de x: os valores esperados para a média e o desvio padrão são 1 e 1
c(mean(x),sd(x))
## [1] 0.9915828 1.0382262
# agora com rexp: esta função gera números aleatórios distribuídos como uma pdf exponencial
y=rexp(1000)
c(mean(y),sd(y))
## [1] 1.009002 0.962707
hist(x,xlim=c(0,5),ylim=c(0,200),breaks=50,col=rgb(1,1,0,0.6),main="distribuição exponencial",xlab="x,y",ylab="frequência")
par(new=TRUE)
hist(y,xlim=c(0,5),ylim=c(0,200),breaks=50,col=rgb(0,1,1,0.3),main="",xlab="",ylab="")

# note que dois bins verdes correspondem a um amarelo

2. integração por Monte Carlo

MC oferece uma forma simples para se integrar uma função positiva \(f(x)\): \[ I = \int_a^b f(x) dx,~~~~f(x)>0\]

Algoritmo do método da rejeição: vamos supor que a função \(f(x)\) possa ser ``coberta’’ por uma função mais simples \(g(x)\) (por exemplo \(g(x)=1\) no intervalo [0,1]), de área \(A\) no intervalo de integração

  • sorteamos N pontos uniformemente dentro de A
  • contamos quantos pontos cairam “dentro” de \(f(x)\), \(N_{ac}\)
  • integral: \[ I \simeq {N_{ac} \over N} A\]
Sobol, A Primer for the Monte Carlo Method (1968…)

Sobol, A Primer for the Monte Carlo Method (1968…)

cálculo de \(\pi\)

Vamos calcular \(\pi\) pela área de um quarto de círculo unitário (\(A=\pi/4\)):

# inicialização
set.seed(123)

# número de pontos na simulação
npt=500

# simulação de x e y
x=runif(npt)
y=runif(npt)

# distância de (x,y) ao centro (0,0)
d=sqrt(x^2+y^2)

# número de pontos 'aceitos'
nac=length(which(d<=1)) 

# valor de pi na simulação: 4 vezes a fração de pontos dentro da área do quarto de círculo
meu_pi=4*nac/npt 
meu_pi
## [1] 3.16
#visualização
plot(x[which(d<=1)],y[which(d<=1)],xlab="x",ylab="y",main=sprintf(paste("pi por Monte Carlo ", npt," pontos")),col='red',asp=1,pch=20)
points(x[which(d>1)],y[which(d>1)],col='blue',pch=20)
xl=seq(0,1,0.01)
yl=sqrt(1-xl^2)
lines(xl,yl,col='green',lwd=2)
abline(v=c(0,1))
abline(h=c(0,1))

Vamos examinar agora a evolução do valor calculado de pi com o número de pontos na simulação:

# função que calcula pi
calc_pi = function(npt){
  c = rep(0,npt)
  Nac = 0
  for(i in 1:npt){
    x = runif(2,0,1)
    if((x[1]*x[1] + x[2]*x[2]) <= 1){
      Nac = Nac + 1
    }
    prop = Nac / i
    meupi = prop *4
    c[i] = meupi
  }
  return(c)
}
npt=1000
res = calc_pi(npt)
ini = 1
plot(res[ini:npt],xlab = 'número de pontos',ylab="pi calc",main="valor calculado de pi vs número de pontos", type = 'l')
lines(rep(pi, npt)[ini:npt], col = 'red')

3. universo homogêneo e isotrópico: distribuição uniforme de pontos em uma esfera

  • objetivo: simular um universo simples: esférico, cartesiano, raio R=1, uniformemente povoado por \(N\) galáxias
  • vamos simular pontos dentro de uma esfera de raio 1, com distribuição angular aleatória usando coordenadas esféricas
  • questão: como simular as coordenadas \((r, \theta, \phi)\) de uma galáxia?
  • regra de ouro: faça um modelo probabilístico que inclua todas as probabilidades relevantes

defina explicitamente todas as probabilidades envolvidas:

parâmetros: N galáxias dentro de um raio R

densidade media: \(\bar n = 3N/(4 \pi R^3)\)

simulação das coordenadas esféricas \((r,\theta, \phi)\):

  • para a galáxia \(i\) geramos 3 números aleatórios uniformemente distribuídos entre 0 e 1: \[ \gamma_{r,i}, \gamma_{\theta,i}, \gamma_{\phi,i}\]

e obtemos valores para \((r_i,\theta_i, \phi_i)\) como: \[r_i = R \gamma_{r,i}^{1/3}\] \[\theta_i = {\rm acos}(1- 2 \gamma_{\theta,i})\] \[\phi_i = 2 \pi \gamma_{\phi,i}\]

# biblioteca gráfica
#install.packages("scatterplot3d")
library(scatterplot3d)

set.seed(1)

# parâmetros do modelo
R=1
N=1000

#densidade media:
dm=3*N/(4*pi*R^3)

# coordenadas esfericas: (r,teta,fi)
teta=acos(1-2*runif(N))
fi=2*pi*runif(N)
r=R*(runif(N))^(1/3)

# coordenadas cartesianas
x=r*sin(teta)*cos(fi)
y=r*sin(teta)*sin(fi)
z=r*cos(teta)

# visualizacao
scatterplot3d(x,y,z,main=" ",asp=1)

Vamos ver se o universo é mesmo uniforme: vamos calcular o número esperado de galáxias dentro de um raio r: \[ N(<r) = \bar n \Big(\frac{r}{R}\Big)^3\]

par(mfcol=c(1,1)) 
#ordeno os raios, do menor para o maior
rs=sort(r)
#densidade media dentro do raio r
nn=seq(1:N)
dens.media=3*nn/(4*pi*rs^3)
plot(rs,nn,xlab="raio r",ylab="número cumulativo de galáxias",type='l',col='blue',main = 'simulado x esperado')
#numero de galaxias esperado dentro do raio r:
rr = seq(0,R,0.01)
nt=N*(rr/R)^3
lines(rr,nt,col="red")

Se voces repetirem algumas vezes a simulação com diferentes sementes iniciais (façam isso com N=100), verão que os gráficos acima variam bastante: isso é um exemplo da chamada “variância cósmica”!

4. MCMC

objetivo do MCMC: produzir amostras de uma distribuição de probabilidades \(P(x)\), onde \(x\) pode ser um vetor

exemplo: MCMC em uma dimensão- amostragem de uma gaussiana de média 2 e desvio padrão 2 (equivalente a rnorm(mu=2,sd=2))

Em geral, quando se trata de probabilidades, é melhor trabalhar com \(\log P(x)\) do que com \(P(x)\) diretamente. Vamos considerar primeiro o caso direto e depois o uso do logaritmo.

Vamos gerar amostras por MCMC para a distribuição de probabilidades \(P(x)=N(x,\mu=2,\sigma=2)\):

O primeiro passo é definir \(P(x)\):

p = function(x){return(dnorm(x, mean = 2, sd = 2, log = FALSE))}

e, em seguida, a função proposta: vamos supor que as propostas sejam também geradas por uma gaussiana, mas de média \(x\) e desvio padrão 1:

funcao_proposta = function(x){return(rnorm(1,mean = x, sd= 1))}

parâmetros do algoritmo de Metropolis:

  • valor inicial para \(x\): \(xini\)
  • número de iterações desejadas: \(niter\)
# função que vai gerar a cadeia via MCMC com o algoritmo de Metropolis-Hastings
mcmc = function(xini, niter){
# defino o vetor que conterá os valores da cadeia  
    cadeia = array(niter+1)
# inicialização da cadeia    
    cadeia[1] = xini
# loop para gerar a cadeia    
    for (i in 1:niter){
# obtenho uma proposta      
        proposta = funcao_proposta(cadeia[i])
# razão entre o valor proposto nesta iteração e o valor anterior da cadeia        
        probab = p(proposta)/p(cadeia[i])
# gero um número aleatório uniformemente distribuído entre 0 e 1 
      runif1=runif(1)
# testo se a proposta é aceita ou não     
        if (probab > runif1) {
# aceita          
            cadeia[i+1] = proposta}
        else{
# não aceita: usa o valor anterior da cadeia         
            cadeia[i+1] = cadeia[i]}
    }
# retorna a cadeia    
    return(cadeia)
}

Esta função devolve uma cadeia com \(niter+1\) valores \(x_i\).

Exemplo:

#inicialização
set.seed(1234)

# vamos começar "longe" do valor médio esperado
xini=-10
niter=1000

#mcmc
y = mcmc(xini,niter)

#visualização da cadeia:
par(mfrow = c(1,2))
x=seq(1,niter+1,1)
plot(x,y,type='l',col='red',ylab='x',xlab='iteração',main='cadeia')

# adiciona uma reta horizontal no valor esperado
abline(h = 2, col="blue" )

hist(y,col='green',main='distribuição de x',xlab='x',ylab='frequência')

# adiciona uma reta vertical no valor esperado
abline(v = 2, col="blue" )

log de P(x)

Muitas vezes é conveniente se trabalhar com o logaritmo da distribuição de probabilidades. Vamos então considerar MCMC com o log de \(P(x)\). Notem a diferença na definição da função:

p = function(x){return(dnorm(x, mean = 2, sd = 2, log = TRUE))}

e no algoritmo mcmc:

mcmc1 = function(xini, niter){
    cadeia = array(niter+1)
    cadeia[1] = xini
    for (i in 1:niter){
        proposta = funcao_proposta(cadeia[i])
# como estou usando o log(p):        
        probab = exp(p(proposta) - p(cadeia[i]))
      runif1=runif(1)
        if (probab > runif1) {
            cadeia[i+1] = proposta}
        else{
            cadeia[i+1] = cadeia[i]}
    }
    return(cadeia)
}

e, rodando de novo:

#inicialização
set.seed(1234)
xini=-10
niter=1000

#mcmc
y = mcmc1(xini,niter)

#visualização da cadeia:
par(mfrow = c(1,2))
x=seq(1,niter+1,1)
plot(x,y,type='l',col='red',ylab='x',xlab='iteração',main='cadeia')
abline(h = 2, col="blue" )
hist(y,col='green',main='distribuição de x',xlab='x',ylab='frequência')
abline(v = 2, col="blue" )

Burn-in

Um certo número de cadeias iniciais (denominadas burn-in) devem ser removidas da análise por não terem ainda atingido uma “região de convergência”. Por exemplo, a figura abaixo mostra as 200 iterações iniciais de 4 cadeias com condições iniciais diferentes:

#inicialização
ncadeia=4
niter=200
xinic=c(-10,0,10,50)

#inicialização de variáveis
y=rep(0,niter+1)
Y=matrix(nrow=niter+1,ncol=ncadeia)

#mcmc
for (i in 1:ncadeia){
y=mcmc1(xinic[i],niter)

# combina  as colunas do vetor y (cadeias)
Y[,i]=cbind(y)
}

#visualização da cadeia:
par(mfrow = c(1,1))
x=seq(1,niter+1,1)
plot(x,Y[,1],type='l',col='red',ylab='x',xlab='iteração',ylim=c(-11,51))
lines(x,Y[,2],type='l',col='green',ylab='x',xlab='iteração')
lines(x,Y[,3],type='l',col='blue',ylab='x',xlab='iteração')
lines(x,Y[,4],type='l',col='black',ylab='x',xlab='iteração')
abline(h = 2, col="blue" )

Nesse exemplo, vemos que as primeiras ~150 iterações devem ser descartadas como burn-in.

Análise

Vamos supor que fizemos uma ou mais simulações e a cadeia convergiu.

Como podemos exibir os resultados? A primeira coisa a se fazer é definir quanto vamos eliminar da cadeia como burn-in. Com o resto da cadeia podemos:

  • determinar a distribuição de \(x\)
  • saber a taxa de aceitação do algoritmo, isto é, a fração de propostas que foram aceitas
  • calcular estatísticas de interesse: média, mediana, variância, quartis…

Vamos fazer uma simulação:

#inicialização
set.seed(1234)
xini=-10
niter=10000

#mcmc
y = mcmc1(xini,niter)

#visualização da cadeia:
par(mfrow = c(1,2))
x=seq(1,niter+1,1)
plot(x,y,type='l',col='red',ylab='x',xlab='iteração',main='cadeia')
abline(h = 2, col="blue" )
hist(y,col='green',main='distribuição de x',xlab='x',ylab='frequência')
abline(v = 2, col="blue" )

Vamos jogar fora, conservadoramente, as 1000 primeiras iterações como burn-in:

burnin = 1000
#aceitacao
aceitacao = 1-mean(duplicated(y[-(1:burnin)]))
aceitacao
## [1] 0.8403511
#media, mediana, desvio padrão
mean(y[-(1:burnin)])
## [1] 2.062565
median(y[-(1:burnin)])
## [1] 2.094431
sd(y[-(1:burnin)])
## [1] 2.002334
#intervalo de confiança de 95%
quantile(y[-(1:burnin)], probs = c(0.025,0.975))
##      2.5%     97.5% 
## -1.983647  5.904849

Convergência

Vamos agora analisar a convergência da cadeia. Para isso vamos jogar fora as iterações de burn-in; em seguida dividimos a cadeia restante em 4 segmentos, calculamos a média e a dispersão dos 4 segmentos e avaliamos.

seg=(niter-burnin)/4
media.seg = c(0,0,0,0)
sd.seg = c(0,0,0,0)
for (i in 1:4){
i1=burnin+1+(i-1)*seg
i2=i1+seg
media.seg[i] = mean(y[i1:i2])
sd.seg[i] = sd(y[i1:i2])
}
#media de cada segmento da cadeia:
media.seg
## [1] 2.109344 2.153425 1.883418 2.105790
#desvio padrão de cada segmento da cadeia
sd.seg
## [1] 1.970430 1.881251 2.054113 2.087099
#media das cadeias
mean(media.seg)
## [1] 2.062994
#desvio padrão das cadeias
sd(media.seg)
## [1] 0.1216623

Um outro procedimento muito comum (principalmente em problemas complexos) é rodar várias cadeias e compará-las, para avaliar a convergência.

Exemplo: vamos comparar resultados de 4 cadeias

#inicialização
ncadeia=4
niter=11000
xinic=c(-1,0,3,5)
#inicialização de variáveis
y=rep(0,niter+1)
Y=matrix(nrow=niter+1,ncol=ncadeia)
#mcmc
for (i in 1:ncadeia){
y=mcmc1(xinic[i],niter)
Y[,i]=cbind(y)}
#visualização da cadeia:
par(mfrow = c(1,1))
x=seq(1,niter+1,1)
plot(x,Y[,1],type='l',col='red',ylab='x',xlab='iteração')
lines(x,Y[,2],type='l',col='green',ylab='x',xlab='iteração')
lines(x,Y[,3],type='l',col='blue',ylab='x',xlab='iteração')
lines(x,Y[,4],type='l',col='black',ylab='x',xlab='iteração')
abline(h = 2, col="blue" )

Podemos, também, calcular algumas estatísticas para essas cadeias (assumindo aqui o mesmo burn-in de antes para todas):

# media de cada cadeia:
media=rep(0,ncadeia)
for (i in 1:ncadeia){media[i]=mean(Y[-(1:burnin),i])}
media
## [1] 1.938203 1.940195 1.981517 2.003466
# desvio padrão de cada cadeia:
dp=rep(0,ncadeia)
for (i in 1:ncadeia){dp[i]=sd(Y[-(1:burnin),i])}
dp
## [1] 2.089776 2.056148 2.138996 1.971478

e, se quisermos, podemos combinar as diferentes cadeias numa “super-cadeia”:

# super cadeia: combinação de todas as cadeias sem o burn-in
Ysuper=c(Y[-(1:burnin),1],Y[-(1:burnin),2],Y[-(1:burnin),3],Y[-(1:burnin),4])
# super media e desvio padrão
mean(Ysuper)
## [1] 1.965845
sd(Ysuper)
## [1] 2.065111

Simulações adicionais

Vamos fazer várias simulações para vocês irem pegando o jeito e verem algumas técnicas de programação usadas em R.

5. Distribuição exponencial e distribuição gama

https://web.stanford.edu/class/bios221/labs/simulation/Lab_3_simulation.html

Pode-se mostrar que a soma de amostras com distribuição exponencial dá uma distribuição gama. Neste exercício vamos gerar 5 amostras com distribuição exponencial com taxa 0.1 e somá-las. Vamos repetir isso 10^5 vezes para ver a cara da distribuição.

# reprodutibilidade
set.seed(123)

# número de exponenciais:
nexps <- 5L
# o sufixo L depois do número significa que o número é inteiro

# taxa da exponencial (valor médio esperado = 1/taxa)
taxa <- 0.1

# qual é o valor esperado da somma dessas exponenciais?

# para gerar 5 amostras com distribuição exponencial e somá-las pode-se usar o comando:
sum(rexp(n=nexps, rate=taxa))
## [1] 28.36911
# vamos repetir isso 10000 vezes
nsim = 10^5L

# para isso vamos usar a função replicate() 
x1 <- replicate(nsim, sum(rexp(n=nexps, rate=taxa))) 

# Para repetições sempre se pode usar loops com o comando 'for', mas isso é muito mais lento que usar funções como replicate() ou apply().

# vamos ver os primeiros valores dessa simulação (comparem com o valor médio esperado)
head(x1)
## [1] 35.31385 23.31460 75.23697 58.10870 53.02865 67.59931
# a distribuição de x1 deve ser igual a função Gama(forma,escala), com forma=nexps e escala=1/taxa
# para fazer uma figura bonita vamos usar a biblioteca ggplot e comparar a distribuição de x1 com a função gama:
require(ggplot2)
## Loading required package: ggplot2
ggplot(data.frame(x1), aes(x1)) + 
  geom_histogram(aes(y=..density..)) +
  stat_function(fun=function(x)dgamma(x, shape=nexps, scale=1/taxa),
                color="red", size=2)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

# essa mesma simulação poderia ser feita de outras maneiras, todas com o mesmo resultado
# por exemplo, com a função sapply()
set.seed(123)
x1 <- sapply(1:nsim, function(i){sum(rexp(n=nexps, rate=taxa))})
head(x1)
## [1] 28.36911 35.31385 23.31460 75.23697 58.10870 53.02865
# ou com a função lapply()
set.seed(123)
x1 <- lapply(1:nsim, function(i){sum(rexp(n=nexps, rate=taxa))})

# quando se aplica uma função simples, como a soma, o jeito mais rápido, em geral, é fazer uma matriz com todas as simulações e então aplicar  a função à natriz. Nesse caso, funções como rowSums() e colSums() são particularmente eficientes.
# usando a função apply()  em uma matriz com nexps linhas e nsim colunas
set.seed(123)
x1 <- apply(matrix(rexp(n=nexps*nsim, rate=taxa), nrow=nexps),2,sum) # o 2 indica que a soma será feita nas colunas
head(x1)
## [1] 28.36911 35.31385 23.31460 75.23697 58.10870 53.02865
# agora usando colSums()
set.seed(123)
x1 <- colSums(matrix(rexp(n=nexps*nsim, rate=taxa), nrow=nexps))

# se seu processador é multicore, voce pode melhorar a velocidade de processamento por paralelização. O pacote parallel tem a função mclapply() que é similar a lapply(). Por default ela use apenas 1 processador, mas você pode mudar mc.cores se seu computador tem mais cores.
require(parallel)
## Loading required package: parallel
# meu laptop:
mc.cores = 4
set.seed(123)
x1 <- mclapply(1:nsim, function(i){sum(rexp(n=nexps, rate=taxa))})

6. Use Monte Carlo para estimar a área entre as curvas da figura

https://www.datasciencecentral.com/profiles/blogs/solving-simple-probability-problems-with-simulation-in-r

x1=seq(0,1,0.1)
y1=x1^2
y2=seq(0,1,0.1)
x2=y2^2
par(mfrow = c(1,1))
plot(x1,y1,xlim=c(0,1),ylim=c(0,1),col='red',lwd=3,type='l',xlab='x',ylab='y')
lines(x2,y2,col='blue',lwd=3)

# reprodutibilidade
set.seed(123)

# número de simulações
n <- 10^6L

# gero n pares (x,y) uniformemente distribuídos entre 0 e 1
df <- data.frame(x = runif(n), y = runif(n)) 

# verifico quais pontos caem entre as curvas
indices <- which(((df$y)^2<=df$x)&((df$x)^2<=df$y)) 

# resultado da simulação
print(paste('Área via  integração por Monte-Carlo = ', length(indices) / n))
## [1] "Área via  integração por Monte-Carlo =  0.333124"
# comparação com a integração numérica:
res <- integrate(function(x) sqrt(x) - x^2,0,1) 
print(paste('Actual Area = ', res$value, ' with absolute error', res$abs.err))
## [1] "Actual Area =  0.333333408167678  with absolute error 7.80933322530332e-05"

7. simulação por MCMC de uma distribuição exponencial

\[P(x) = exp(-x)\]

# reprodutibilidade
set.seed(1234)

# função a ser amostrada
p = function(x){return(dexp(x, log = FALSE))}

# a função proposta deve propor valôres dentro do espaço amostral da distribuição de probabilidades considerada. A distribuição exponencial é definida para x> 0, enquanto que uma exponencial pode ter valores negativos. Assim, precisamos modificar a função proposta gaussiana para ela propor apenas valores positivos:
funcao_proposta = function(x){
a = -1
while(a < 0){a = rnorm(1,mean = x, sd= 1)}
return(a)}

mcmc = function(xini, niter){
# defino o vetor que conterá os valores da cadeia  
    cadeia = array(niter+1)
# inicialização da cadeia    
    cadeia[1] = xini
# loop para gerar a cadeia    
    for (i in 1:niter){
# obtenho uma proposta      
        proposta = funcao_proposta(cadeia[i])
# razão entre o valor proposto nesta iteração e o valor anterior da cadeia        
        probab = p(proposta)/p(cadeia[i])
# gero um número aleatório uniformemente distribuído entre 0 e 1 
      runif1=runif(1)
# testo se a proposta é aceita ou não     
        if (probab > runif1) {
# aceita          
            cadeia[i+1] = proposta}
        else{
# não aceita: usa o valor anterior da cadeia         
            cadeia[i+1] = cadeia[i]}
    }
# retorna a cadeia    
    return(cadeia)
}

# vamos começar "longe" do valor esperado
xini=10
niter=2000

#mcmc
y = mcmc(xini,niter)

#visualização da cadeia:
par(mfrow = c(1,2))
x=seq(1,niter+1,1)
plot(x,y,type='l',col='red',ylab='x',xlab='iteração',main='cadeia')
# adiciona uma reta horizontal
abline(h = 1, col="blue" )
hist(y,col='green',main='distribuição de x',xlab='x',ylab='frequência',breaks=20)
# adiciona uma reta vertical
abline(v = 1, col="blue" )

8. random walk em duas dimensões

Vamos considerar um modelo de passeio aleatório (random walk) discreto no plano (x,y).

  • Vamos supor que cobrimos o plano (x,y) com uma malha de separação unitária
  • Vamos supor que, em um dado ponto (x,y) da malha , só 4 movimentos (“passos”) são possíveis: dx=1 dx=-1 dy=1 dy=-1
  • Começamos na origem e a partir daí escolhemos movimentos aleatórios

Qual é a distância média à origem depois de 1000 passos? Faça 5000 simulações para calcular esta distância média.

Pode-se mostrar que a distribuição esperada das distâncias finais, no caso de de N passos (\(N >> 1\)) é dada pela distribuição de Rayleigh: \[ P(d|N)=\frac{2d}{N}e^{-d^2/N} \] e que seu valor esperado é \[ \bar d = \sqrt{N} \]

# reprodutibilidade
set.seed(1234)

#parametros da simulacao
nsim=5000
npassos=1000

# definindo o vetor que conterá a distância final da nave em cada simulação
d=rep(0,nsim)
for (isim in 1:nsim) {
x=rep(0,npassos)
y=rep(0,npassos)
for (n in 2:npassos) {
x[n]=x[n-1]
y[n]=y[n-1]
r=runif(4)  # gero 4 numeros para escolher o passo
imax=which.max(r)  # passo mais provável
if (imax==1) {x[n]=x[n]+1}
if (imax==2) {x[n]=x[n]-1}
if (imax==3) {y[n]=y[n]+1}
if (imax==4) {y[n]=y[n]-1}
}
# distância final nesta simulação:
d[isim]=(x[npassos]^2+y[npassos]^2)^0.5
}

plot(x,y,type="l",main="exemplo de uma trajetória",col="red")
# ponto de partida
points(0,0,col='blue',pch=20)

hist(d,main="histograma das distâncias finais",col='salmon')

# distância média:
mean(d)
## [1] 27.98005
summary(d)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00   17.12   26.40   27.98   37.22   99.04
#valor esperado de d: npassos^1/2
sqrt(npassos)
## [1] 31.62278
# tempo de processamento:
Sys.time()-t00
## Time difference of 16.24126 secs