Prática opcional - Simulando o passeio do bêbado
Imagine que um indivíduo sai do Natividade bastante embriagado. Ele bebeu tanto que está “trocando as pernas” na rua, dando passos para a direita e para a esquerda com igual probabilidade. Em estatística, um fenômeno assim é conhecido como “random walk” ou passeio aleatório e suas aplicações são inúmeras em várias áreas da ciência (além de prever quando nosso amigo vai cair na sarjeta). Para simular um random walk, precisamos de alguns componentes.
Teste o seguinte comando no console do R
rbinom(1, 1, 0.5)## [1] 0Usamos esse comando na aula passada na simulação de jogo de “cara ou coroa”. No nosso caso, esse comando gera como resposta o número 0 ou o número 1 com a mesma probabilidade. No enunciado acima, nosso amigo dá passos à direita ou à esquerda com a mesma probabilidade. Vamos usar o comando acima em algo um pouco mais complexo. Imagine que onde nosso amigo está na rua é descrito por um objeto chamado posicao e que ele começa no meio da rua, um ponto que arbitrariamente chamamos de zero (ou seja, a posicao[1] é 0) e dá um passo à direita (posicao[1] é 1) ou esquerda (posicao[1] é -1) com a mesma probabilidade.
posicao <- numeric()
posicao[1] <- 0
posicao[2] <- ifelse(rbinom(1, 1, 0.5) == 0, posicao[2] <- posicao[1] + 1, posicao[2] <- posicao[1] -
1)Há várias novidades aqui. A primeira é que temos um comando em várias linhas. Poderíamos na verdade escrever tudo em uma linha só, mas como você pode ver, fica muito mais simples para entender se separamos em linhas diferentes. Uma alternativa que facilitará muito a sua vida é abrir uma novo arquivo no R (arquivo>novo ou algo parecido), escrever todos os comandos e depois enviá-los todos de uma vez ao console. Em aula eu ensinarei como fazer isso. Alternativamente, você pode simplesmente selecionar tudo, copiar e colar no console.
Agora vamos olhar os comandos passo a passo. O primeiro é necessário para dizer ao R para criar um objeto numérico. Não podemos querer atualizar a posição do nosso amigo se não ensinarmos ao R o que queremos dizer com posição. No segundo comando, ocupamos a primeira “prateleira” do objeto posicao com a condição inicial - nesse caso, a posição 0.
O terceiro comando é um pouco mais complexo. Usamos o comando ifelse. Este comando funciona da seguinte maneira: ifelse(algum teste que eu queira, o que fazer se o teste for verdadeiro, o que fazer se o teste for falso). No nosso caso, usamos o primeiro comando que aprendemos hoje para gerar um número aleatório (0 ou 1) e perguntamos ao R se o número simulado é zero. Se sim, o R faz o que está na primeira linha abaixo. Se não, ele faz o que está na segunda linha abaixo. Ou seja, se o número simulado é maior que 0.5, o nosso amigo dá um passo à direita (sua posição 2 é 0+1=1) ou à esquerda (sua posição 2 é 0-1=-1).
Como poderíamos simular este passeio depois de 5 passos? A maneira mais trivial seria repetir o comando ifelse 5 vezes.
posicao <- numeric()
posicao[1] <- 0
posicao[2] <- ifelse(rbinom(1, 1, 0.5) == 0, posicao[2] <- posicao[1] + 1, posicao[2] <- posicao[1] -
1)
posicao[3] <- ifelse(rbinom(1, 1, 0.5) == 0, posicao[3] <- posicao[2] + 1, posicao[3] <- posicao[2] -
1)
posicao[4] <- ifelse(rbinom(1, 1, 0.5) == 0, posicao[4] <- posicao[3] + 1, posicao[4] <- posicao[3] -
1)
posicao[5] <- ifelse(rbinom(1, 1, 0.5) == 0, posicao[5] <- posicao[4] + 1, posicao[5] <- posicao[4] -
1)
posicao[6] <- ifelse(rbinom(1, 1, 0.5) == 0, posicao[6] <- posicao[5] + 1, posicao[6] <- posicao[5] -
1)
posicao## [1] 0 1 0 1 2 1Você pode visualizar este passeio também.
plot(posicao, type = "l")
plot of chunk unnamed-chunk-4
E se quisermos saber como esse passeio ficaria depois de 1000 passos? Escrever passo a passo não parece muito eficiente. Aprenderemos outra ferramenta muito útil no R: os loops. Primeiro vamos fazer um teste para aprendermos como o loop funciona:
test <- numeric()
test[1] <- 13
for (i in 2:10) {
test[i] <- test[i - 1] + 1
}O que esse comando fez? Ele usou o valor em [1] para calcular o valor de [2], depois usou o valor de [2] para calcular o valor de [3], etc. O loop tem a seguinte estrutura. Primeiro temos (i in 2:10). Isso quer dizer que o R vai repetir os comandos substituindo i por 2, 3, 4, 5, 6, 7, 8, 9 e 10 no comando que está entre {}. Ou seja, ele repete o comando que está entre {} e o repete, substituindo onde está i pelo respectivo número. Voltando ao passeio aleatório, podemos ter:
posicao <- numeric()
posicao[1] <- 0
for (i in 2:1000) {
posicao[i] <- ifelse(runif(1, min = 0, max = 1) > 0.5, posicao[i] <- posicao[i -
1] + 1, posicao[i] <- posicao[i - 1] - 1)
}Podemos visualizar o passeio depois de 1000 passos:
plot(posicao, type = "l")
plot of chunk unnamed-chunk-7
Agora digamos que você gostou tanto dos comandos que você fez acima que você quer repetir várias vezes. Mas é desajeitado ter que repetir tudo cada vez. Novamente o R tem uma solução: faremos nossa própria função!
passeio <- function(passos) {
posicao <- numeric()
posicao[1] <- 0
for (i in 2:passos) {
posicao[i] <- ifelse(rbinom(1, 1, 0.5) == 0, posicao[i] <- posicao[i -
1] + 1, posicao[i] <- posicao[i - 1] - 1)
}
return(posicao)
}Parabéns! Agora que você ensinou ao R o que fazer quando pedir pra ele fazer passeio, basta fazer:
passeio(100)## [1] 0 1 0 1 0 1 0 -1 -2 -3 -4 -3 -4 -3 -2 -3 -2
## [18] -3 -2 -3 -4 -5 -4 -5 -6 -5 -4 -5 -6 -7 -6 -7 -6 -7
## [35] -8 -7 -8 -9 -10 -9 -8 -9 -10 -9 -10 -9 -8 -7 -8 -9 -10
## [52] -11 -10 -9 -8 -9 -8 -9 -8 -9 -8 -7 -6 -7 -6 -5 -4 -5
## [69] -4 -5 -4 -3 -2 -3 -4 -5 -6 -7 -8 -9 -8 -9 -10 -11 -10
## [86] -9 -10 -9 -10 -11 -10 -11 -10 -9 -8 -9 -10 -11 -12 -11Agora que estamos empolgados, que tal simularmos 10 caminhantes ao mesmo tempo, com cada um dando mil passos? Como já fizemos a função passeio, isso é muito fácil:
resultados <- replicate(10, passeio(1000))Agora “abra” o objeto resultados. Você verá uma matriz onde cada coluna é um caminhante. Note que uma matriz é um tipo de objeto diferente dos que havíamos encontrado. Agora há linhas e colunas. Por exemplo, você pode se direcionar ao valor na linha 2, coluna 3 usando a notação objeto[2,3]. Se eu quero todos os valores na primeira coluna, basta deixar vazio o campo das linhas, como objeto[,1].
Podemos plotar tudo ao mesmo tempo.
plot(resultados[, 1], type = "l", ylim = c(min(resultados), max(resultados)))
lines(resultados[, 2])
lines(resultados[, 3])
lines(resultados[, 4])
lines(resultados[, 5])
lines(resultados[, 6])
lines(resultados[, 7])
lines(resultados[, 8])
lines(resultados[, 9])
lines(resultados[, 10])
lines(rowSums(resultados)/10, col = "red")
plot of chunk unnamed-chunk-11
Há diversos fenômenos biológicos que podem ser modelados com a abordagem acima, como a deriva genética em populações naturais.
As simulações que realizamos são realizadas essencialmente em duas dimensões, já que nosso amigo pode caminhar somente para a direita ou esquerda. Por outro lado, podemos imaginar exemplos onde passeios aleatórios em duas dimensões sejam de interesse. Por exemplo, podemos imaginar um conjunto de indivíduos de uma espécie invasora que se originam em um ponto de entrada em um país e podemos estar interessados em como estes indivíduos se disseminam geograficamente. Para isso precisaremos de coordenadas X e Y. Para facilitar nossa vida, imaginemos um mundo quadriculado onde em cada momento um indivíduo escolhe aleatoriamente se vai para o norte/sul ou leste/oeste. Como o deslocamento horizontal e vertical são independentes, podemos utilizar a mesma função acima duas vezes:
X <- passeio(10000)
Y <- passeio(10000)
plot(X, Y, type = "l")
plot of chunk unnamed-chunk-12
Agora digamos que sabemos que a população inicial era de 1000 indivíduos e que se movimentam uma vez por dia. Como poderíamos esperar que estes indivíduos se distribuiriam no espaço desde sua introdução até um ano depois?
simX <- replicate(1000, passeio(365))
simY <- replicate(1000, passeio(365))Isso deve demorar alguns minutos. Agora podemos plotar esses dados ao longo do tempo. Pela maneira que o R plota gráficos, é mais simples começar pelo final:
plot(simX[360, ], simY[360, ], xlab = "X", ylab = "Y", col = "yellow") # posições após quase 1 ano
points(simX[270, ], simY[270, ], xlab = "X", ylab = "Y", col = "orange") # posições após 9 meses
points(simX[180, ], simY[180, ], xlab = "X", ylab = "Y", col = "red") # posições após 6 meses
points(simX[90, ], simY[90, ], xlab = "X", ylab = "Y", col = "darkred") # posições após 3 mesesExercício avançado: * Repita o terceiro gráfico utilizando um loop no lugar de plotar individualmente cada série.