Introdução

Neste documento faremos algumas simulações e análises dos lançamentos de Moeda.

Vamos começar carregando os pacotes necessários.

#Descomente se for necessário instalar os pacotes a seguir:
#install.packages("dplyr")
suppressMessages(library(dplyr))

Lançando Moedas

Vamos começar criando uma função flips. A função flips(v) simula o lançamento de v moedas. As saídas de nossas moedas serão \(\{0,1\}\).

sample.space <- c(-1,1)
theta <- 0.5 # moeda honesta
flips <- function(v){ sample(sample.space, 
                size = v, 
                replace = TRUE, 
                prob = c(theta, 1 - theta))}
# Jogar uma moeda 100 vezes
flips(100)
##   [1]  1 -1 -1  1 -1  1  1  1 -1  1  1  1 -1  1  1  1  1  1  1 -1 -1  1  1
##  [24]  1 -1 -1 -1 -1 -1  1 -1 -1  1 -1 -1  1 -1  1  1 -1 -1 -1 -1  1  1 -1
##  [47]  1  1 -1  1  1  1 -1 -1 -1  1  1 -1  1  1  1 -1  1 -1  1 -1  1  1  1
##  [70] -1  1  1 -1  1  1 -1 -1 -1  1 -1  1  1 -1 -1  1  1 -1 -1  1  1  1 -1
##  [93]  1 -1  1  1 -1 -1 -1 -1

Vejamos quantas caras e quantas coroas aparecem em 10000 lancamentos

table(flips(10000))
## 
##   -1    1 
## 4913 5087

Somas Cumulativas e Passeios Aleatórios

Vamos calcular as somas parciais das primeiras \(k\) entradas para \(k=1,\dots n\). Ou seja se denotarmos por \(X_i\) o resultado do \(i\)-ésimo lançamento. A soma até \(k\) é definida como \[ S_k=\sum_{i=1}^k X_i\]

A soma \(S_k\) é o que chamamos de passeio aleatório unidimensional.

Vejamos as somas parciais.Isso pode ser feito com a função cumsum

temp<-flips(30)
temp
##  [1]  1 -1  1 -1 -1 -1 -1  1  1 -1 -1 -1 -1 -1 -1 -1 -1 -1  1  1  1 -1  1
## [24] -1 -1  1  1  1 -1 -1
cumsum(temp)
##  [1]   1   0   1   0  -1  -2  -3  -2  -1  -2  -3  -4  -5  -6  -7  -8  -9
## [18] -10  -9  -8  -7  -8  -7  -8  -9  -8  -7  -6  -7  -8

Vamos plotar as somas parciais ou cumulativas.

 plot(cumsum(flips(100000)), type = 'l', main="Passeio Aleatório",ylab ="Somas Parciais" )

Vamos simular o comportamento de cinco somas cumulativas no mesmo gráfico

plot(cumsum(flips(100000)), type = 'l',col="red", main="Passeios Aleatórios",ylab ="Somas Parciais",ylim=c(-400,400) )
lines(cumsum(flips(100000)), type = 'l',col="green")
lines(cumsum(flips(100000)), type = 'l',col="blue")
lines(cumsum(flips(100000)), type = 'l',col="pink")
lines(cumsum(flips(100000)), type = 'l',col="yellow")

Teorema Central do Limite

Vamos criar um histograma das posições após 1000 lançamentos para 10000 passeios

 # Criando uma lista vazia e armazenando as somas cumulativas
 resultados <- list()
 for(i in 1:12000) {
     coinTosses   <- cumsum(sample(c(-1,1), 1000, replace = TRUE)) 
     resultados[[i]] <- coinTosses[length(coinTosses)]
 }
 
 # Criar um histograma. Defina um título e defina a cor e quebra
   hist(unlist(resultados), main = "Histograma de todas as posições finais",col = "lightblue", breaks = 100)
 
 #Colocar uma linha vertical em 0 com uma largura de 2 para mostrar a média da distribuição
 abline(v = 0, col = "red", lwd = 2)

Esse gráfico ilustra bem o TCL:

Teorema do Limite Central de Lindeberg–Lévy

Dados \(\{X_1, X_2, …\}\) uma sequência de variáveis aleatórias i.i.d com \(E[X_i] = \mu\) e \(Var[X_i] = \sigma^2 <\infty\). Então o limite quando \(n\) tende ao infinito de \(\frac{S_n-n\mu}{\sigma\sqrt{n}}\) converge em distribuição para a distribuição normal \(N(0,1)\).

\[ \text{ CLT: } \quad \lim_{n\rightarrow \infty} P\Bigl(\frac{S_n-n\mu}{\sigma\sqrt{n}} < x \Bigr) = \Phi(x)\]

Números Normais de Borel

Vamos analisar a média até a posição \(n\). Veja que a média converge rapidamente para 0.

 plot(cummean(flips(1000)), type = 'l', main="Médias Parciais ")

Vamos plotar as médias parciais de \(1,\dots,2000\)

mediaseq<-function(v){mean(flips(v))}
  s<-lapply(1:2000,mediaseq)
  boxplot(s)

Vamos estimar qual a probabilidade que uma sequência possua média maior que \(p\)

mediasfinais <- list()
 for(i in 1:10000) {
    mediasfinais[i]<-mean(flips(1000))
 }
mediasfinais<-data.frame(mediasfinais)

contamaioresque <-function(p){ length(mediasfinais[ ,which(mediasfinais[1,]>p) ])/10000}

# Veja qual a percentagem das médias é maior que 0.1
contamaioresque(0.1)
## [1] 9e-04
maiores<-lapply(1:20/100,contamaioresque)
barplot(t(maiores))

mediaapos1000<-function(v){mean(flips(1000))}
  s<-lapply(1:2000,mediaapos1000)
barplot(t(s), main = "Média após 1000 lançamentos")

Lei Fraca

A lei fraca de grandes números (também chamada de lei de Khintchine) afirma que a média da amostra converge em probabilidade para o valor esperado

\[{\bar {X}}_{n}= S_n/n \xrightarrow{P} \mu \qquad \text{ quando }\ n\to \infty . \]

#setting a parameters of Bi(n, p)
n <- 1000;
p <- 0.4;
 
#dataframe
df <- data.frame(bi = rbinom(n, 1, p)  ,count = 0, mean = 0);
ifelse(df$bi[1] == 1, df[1, 2:3] <- 1, 0);
## [1] 0
for (i in 2 : n){
  df$count[i] <- ifelse(df$bi[i] == 1, df$count[i]<-df$count[i - 1]+1, df$count[i - 1]);
  df$mean[i] <- df$count[i] / i;
};
 
#graph
plot(df$mean, type='l',
      main = "Simulação da Lei Fraca dos Grandes Números",
      xlab="Numbers", ylab="Sample mean")
abline(h = p, col="red")

Lei Forte dos Grandes Números para Experimentos de Bernoulli

A lei forte de grandes números afirma que a média da amostra converge quase com certeza ao valor esperado

A lei fraca afirma que, dado \(n\) um número suficientemente grande, a média \({\displaystyle {\overline {X}} _ {n}}\) é provável que esteja perto de \(\mu\). Assim, deixa aberta a possibilidade de que \({\displaystyle | {\overline {X}} _ {n} - \mu |> \varepsilon}\) ocorra um número infinito de vezes, embora em intervalos infreqüentes.

A lei forte mostra que isso quase certamente não ocorrerá. Em particular, isso implica que com a probabilidade 1, temos isso para qualquer \(\varepsilon> 0\) a desigualdade \({\displaystyle | {\overline {X}} _ {n} - \mu | <\varepsilon}\) é válido para \(n\) suficientemente grande.

\[{\displaystyle {\begin{matrix}{}\\{\bar {X}}_{n}\ {\xrightarrow {\text{q.c.}}}\ \mu \qquad {\textrm {quando}}\ n\to \infty .\\{}\end{matrix}}}\]

Isto é

\[\displaystyle P \!\left(\lim _{n\to \infty }{\bar {X}}_{n}=\mu \right)=1.\]

par(mfrow=c(3,3),oma=c(0,0,5,0))
layout(matrix(1:9,nrow=3))
sllnBern <- function(n,p) {
for (i in 1:9) {
    seq <- rbinom(n,1,p)
    avgs <- cumsum(seq)/(1:n)
    plot(avgs,type="l",xlab="n",ylab="Média")
    abline(h = p)
    }
 mtext(paste("Lei Forte dos Grandes Números ( n = ",n,")"),side=3, font=2,outer=T)
 mtext(paste("Experimentos de Bernoulli ( p = ",p,")"),side=3, outer=T, font=2, line =-2)
 }
sllnBern(100,.5)

sllnBern(1000,.5)

sllnBern(10000,.5)

Lei dos Logaritmos Iterados

#Definir a função limsup
limsup <- function(x){
 z = rep(0,length(x))
for(i in 1:length(x)){
 z[i] = max(x[i:length(x)])
 }
 return(z)
}
Vamos plotar as somas parciais de $X_1+\dotsX_n$ e analisar seu crescimento
as<- suppressWarnings(function(n){sqrt(2*n*log(log(n)))})
asx<- suppressWarnings(function(n){-sqrt(2*n*log(log(n)))})

mediaseqit<-function(v){sum(flips(v))}
  sit<-as.vector(sapply(seq(5,5000),mediaseqit))
    plot(sit, type = 'l',col="darkgreen", main="Passeios Aleatórios",ylab ="Somas Parciais",ylim=c(-400,400) )
   suppressWarnings(curve(as, add=TRUE,col="red"))
  suppressWarnings(curve(asx, add=TRUE,col="blue"))
  legend("topleft", inset=.025, title="Funções",
   c("2*n*log(log(n)","-2*n*log(log(n)"), fill=c("red","blue"))