• Durante o curso buscamos por padrões nos dados para entender comportamentos ou prever resultados de ações. Nem sempre esses padrões são fáceis de identificar, choques aleatórios costumam afetar as variáveis de interesse e “contaminar” os padrões que observamos.
• Um dos trabalhos de economistas aplicados é encontrar padrões que podem estar ocultos por conta de choques aleatórios, ou seja, de alguma forma temos de separar o que é ruído aleatório e o que é um padrão sistemático dos dados.
• A teoria da probabilidade fornece um conjunto de ferramentas matemáticas que nos permite medir os padrões aleatórios encontrados nos dados.

• Antes de comerçarmos a falar de teoria é importante refletir um pouco sobre o que é probabilidade. O que queremos dizer quando afirmamos que a probabilidade de algo acontecer é 30% ou 70%.
• Existem duas formas de pensar probabilidade: frequentista (objetiva) e bayesiana (subjetiva).

• Qual a probabilidade do Botafogo perder o próximo jogo contra o Flamengo?
  • Considerando que o hino do Botafogo deixa claro que o Glorioso não pode perder para ninguém acredito que a probabilidade é zero. (resposta subjetiva)
  • Procurei os dados na internet, encontrei que em 380 jogos contra o Flamengo o Botafogo ganhou 114 (30%), empatou 124 (32,6%) e perdeu 142 (37,4%). Acredito que a resposta é 37,4%. (resposta objetiva)

• Na estatística tradicional, chamada de frequentista, o conceito de probabilidade está associado à resposta que chamei de objetiva.
• Na estatística bayesiana o conceito de probabilidade está associado à resposta que chamei de subjetiva.

• Apesar das diferenças conceituais, estatísticos frequentistas e bayesianos usam a teoria da probabilidade desenvolvida pelo matemático russo Andrey Kolmogorov no começo do século XX. É essa abordagem matemática, comum a frequentistas e bayesianos, que vamos estudar nesta seção.

• Existem duas formas predominantes de interpretar probabilidade. De acordo com a abordagem frequentista, probabilidade representa o limite da frequência relativa que evento de interesse ocorre quando o número de experimentos repetidos nas mesmas condições vai para infinito. A abordagem bayesiana interpreta probabilidade como uma crença subjetiva a respeito da probabilidade de um evento ocorrer.

• Para definir probabilidade usaremos três conceitos: experimento, espaço amostral e evento.
  • experimento: uma ação ou conjunto de ações que produzem os eventos estocásticos de interesse.
  • espaço amostral: o conjunto de todos os possíveis resultados de um experimento, tipicamente denotado por \(\Omega\).
  • evento: um subconjunto do espaço amostral.

• Exemplo: Cara ou coroa
  • experimento: jogar uma moeda
  • espaço amostral: {cara, coroa}
  • eventos: {cara}, {coroa}

• Exemplo: Botafogo vs Flamengo
  • experimento: jogo disputado pelas duas equipes
  • espaço amostral: {Botafogo ganha, Flamengo ganha, Empate}
  • eventos: {Botafogo ganha}, {Flamengo ganha}, {Empate}

• Pela definição de eventos, subconjunto do espaço amostral, existem vários possíveis eventos, de fato o próprio espaço amostral é um evento.
• Dizemos que um evento ocorreu quando o conjunto que define o evento é o resultado do experimento.

• Considere a decisão de cada eleitor nas eleições de 2022 como um experimento. Para simplifcar considere quatro possiveis ações do eleitor: votar no Lula, votar no Bolsonaro, abstenção (incluindo branco, nulo e não votar) e votar em outro candidato.
• O espaço amostral será: \(\Omega = \{\mbox{votar no Lula}, \mbox{votar no Bolsonaro}, \mbox{abstenção}, \mbox{votar em outro}\}\)

• Podemos considerar vários eventos nesse espaço amostral, por exemplo:
  • {votar no LUla}
  • {votar no Lula, votar no Bolsonaro, votar em outro}, ou seja, não se abster.
  • {votar no Lula, abstenção, votar em outro}, ou seja, não votar no Bolsonaro.
  • {abstenção, votar em outro}, ou seja, não votar nem em Lula nem em Bolsonaro.

• Vamos começar definindo probabilidade para o caso mais simples de todos: cada resultado é igualmente provável de acontecer. Nesse caso, a probabilidade de um evento A, denotada por P(A), pode ser calculada como o número de elementos no conjunto A dividido pelo número de elementos no espaço amostral. \[P(A) = \frac{\mbox{número de elementos de A}}{\mbox{número de elemento de }\Omega}\]

• Como exemplo, considere o experimento de jogar um moeda três vezes. Denote dar cara como H, de heads, e dar coroa como T, de tail. O espaço amostral é: \[\Omega = \{HHH, HHT, HTH, THH, HTT, THT, TTH, TTT\}\]
• Podemos calcular a probabilidade do evento dar pelo menos duas caras contando o número de elementos no conjunto relevante: \[A = \{HHH, HHT, HTH, THH\}\]
• A probabilidade de dar pelo menos duas caras é: \[P(A)=\frac{4}{8}=0,\!5\]

• Toda a moderna teoria de probabilidade é construída a partir de três axiomas, ou seja, a partir desses axiomas pdoemos provar todos os teoremas e regras sobre probabilidade.
• Parece Microeconomia!

• Os axiomas de probabilidade são dados por:
  • A probabilidade de qualquer evento A é não-negativa: \(P(A) \geq 0\)
  • A probabilidade de que um dos resultados no espaço amostral ocorra é um: \(P(\Omega)=1\)
  • Regra da adição, se os eventos A e B são mutuamente exclusivos então \(P(A\mbox{ ou }B)=P(A)+P(B)\)

• Dois eventos A e B são mutuamente exclusivos quando não compartilham de um resultado.
• Os eventos Botafogo ganha e Flamengo ganha são mutuamente exclusivos, desta forma: \[P(\mbox{Botafogo ganhar ou Flamengo ganhar}) = \\P(\mbox{Botafogo ganhar}) + P(\mbox{Flamengo ganhar})\].

• Os eventos Botafogo não ganhar e Flamengo não ganhar não são mutuamente exclusivos pois compartilham o resultado Empate. Note que Botafogo não ganhar = {Flamengo ganhar, Empate} e Flamengo não ganhar = {Botafogo ganhar, Empate}.

• Se A e B não são eventos mutuamente exclusivos vale que: \(P(A\mbox{ ou }B)=P(A)+P(B)-P(A\mbox{ e }B)\).
• No exemplo, P(Botafogo não ganhar ou Flamengo não ganhar) = P(Botafogo não ganhar) + P(Flamengo não ganhar) - P(Empate).
• Repare que a regra acima vale para quaisquer dois eventos, pois se A e B são mutuamente exclusivos então \(P(A\mbox{ e }B)=0\).

• Para quaisquer eventos A e B, a regra da adição é dada por: \(P(A\mbox{ ou }B)=P(A)+P(B)-P(A\mbox{ e }B)\).

• Para um dado evento \(A\) o conjunto dos elementos do espaço amostral que não estão em A é chamado de complemento de A (ou complementar de A) e denotado por \(A^c\).
• Pelo segundo axioma temos que \(P(A)+P(A^c)=1\), que implica em \(P(A^c)=1-P(A)\).
• Por fim, repare que um evento A pode ser decomposto em dois eventos mutuamente exclusivos. A parte de A que é comum a B, \(A\) e \(B\), e a parte de A que não é comum a B, \(A\) e \(B^c\).

• Para quaisquer eventos A e B, a lei da probabilidade total é dada por: \(P(A)=P(A\mbox{ e }B) + P(A\mbox{ e }B^c)\).

• Para calcular a probabilidade do evento A precisamos contar o número de elementos no evento A e o número de elementos no espaço amostral, \(\Omega\).
• No exemplo do “cara ou coroa” contamos o número de elementos fazendo permutações, que faz referência ao número de maneiras que objetos podem arrumados.

• Por exemplo, considere três objetos A, B e C. Existem seis maneiras de arrumar esses objetos: \(\{ABC, ACB, BAC, BCA, CAB, CBA\}\). Como podemos calcular o número de permutações possíveis sem termos de listar todas as combinações?
• Tem um jeito! Existem três escolhas para o primeiro objeto: A, B ou C. Depois que o primeiro objeto for escolhido existem duas escolhas para o segundo objeto. Para o terceiro objeto só existe uma escolha.

• Para calcular o número final de combinações só precisamos multiplicar o número de opções (setas) em cada etapa, ou seja, o número de combinações é \(3\times 2 \times 1=6\).

• O número de permutações quando arrumando \(k\) objetos de um total de \(n\) objetos distintos é denotado por \(_nP_k\) e dado por: \[_nP_k=n\times(n-1)\times\cdots\times(n-k+2)\times(n-k+1)=\frac{n!}{(n-k)!}\]

• No exemplo anterior tínhamos \(n=3\) e \(k=3\), logo o número de permutas é dado por \(_3P_3=\frac{3!}{0!}=6\).
• Outro exemplo é arranjar 4 cartas de um baralho com 13 cartas únicas, a reposta é \(_{13}P_{4}=\frac{13!}{(13-4)}=13\times 12 \times 11 \times 10 = 17160\)

• O problema dos aniversários é um exemplo conhecido e contraintuitivo de permutações. Quantas pessoas devem ter em uma sala para que a probabilidade de pelo menos duas pessoas terem a mesma data de aniversário seja maior do que 50%?
• A resposta é 23, para muita gente parece pouco.

• Para resolver esse problema usando permutação, primeiro repare que: \[P(\mbox{pelo menos duas pessoas tem aniversário no mesmo dia})=\\1-P({\mbox{ninguém tem aniversário no mesmo dia}})\]
• Faz sentido? Pense nos complementos.

• Desta forma podemos responder a questão calculando a probabilidade de ninguém ter nascido no mesmo dia.
• Seja \(k\) o número de pessoas na sala. Para calcular a probabilidade de ninguém ter o mesmo aniversário precisamos contar o número de maneiras que \(k\) pesoas podem ter aniversários diferentes e dividir pelo número de maneiras de selecionar \(k\) possíveis datas de aniversário.

• Desde que toda dada de aniversário seja igualmente provável, o numerador consiste em quantas maneiras podemos ter \(k\) aniversários únicos dentre os 365 dias possíveis. Esse valor é dado por \(_{365}P_k = \frac{365!}{(365-k)!}\)
• Para o denominador considere que cada pessoa pode ter um de 365 dias para aniversário, como são \(k\) pessoas o total de aniversários possíveis é \(365^k\).

• Probabilidade de ninguém ter nascido no mesmo dia é: \[P({\mbox{ninguém nasceu no mesmo dia}})\\=\frac{\mbox{número de maneiras de k dias únicos de aniversário}}{\mbox{número possível de aniversários não únicos para k pessoas}}\\=\frac{\frac{365!}{(365-k)!}}{365^k} = \frac{365!}{365^k(365-k)!}\]

• Desta forma, \[P(\mbox{pelo menos duas pessoas tem aniversário no mesmo dia})\\=1-\frac{365!}{365^k(365-k)!}\]

• Calcular o valor acima não é trivial nem para valores pequenos de \(k\), em casos assim é comum aplicar o logaritmo na expressão. \[\log P({\mbox{ninguém nasceu no mesmo dia}})=\\\log 365! - k\log365 - \log(365-k)!\]
• Depois aplicamos a exponencial em \(\log P({\mbox{ninguém nasceu no mesmo dia}})\) para encontrarmos o valor da probabilidade.

• No R usamos a função factorial() para calcular o fatorial de lfactorial() para calcular o log do fatorial.

factorial(5)

## [1] 120

factorial(10)

## [1] 3628800

factorial(30)

## [1] 2.652529e+32

lfactorial(5)

## [1] 4.787492

log(120)

## [1] 4.787492

exp(lfactorial(5))

## [1] 120

factorial(365)

## [1] Inf

lfactorial(365)

## [1] 1792.332

• Para encontrar a probabilidade de pelo menos duas pessoas terem a mesma data de aniversário em uma turma de \(k\) vamos criar uma função que terá \(k\) como argumento.

library(tidyverse)

birthday <- function(k) {
  logdenom <- k * log(365) + lfactorial(365-k)
  lognumer <- lfactorial(365)
  pr <- 1 - exp(lognumer - logdenom)
  pr
}

• Agora podemos criar um tibble com as probabilidades para diversos tamanhos de turmas.

bday <- tibble(k=1:50, pr=birthday(k))
head(bday)

## # A tibble: 6 × 2
##       k       pr
##   <int>    <dbl>
## 1     1 2.27e-13
## 2     2 2.74e- 3
## 3     3 8.20e- 3
## 4     4 1.64e- 2
## 5     5 2.71e- 2
## 6     6 4.05e- 2

• Por fim, podemos ilustrar o resultado em um gráfico.

ggplot(bday, aes(k, pr)) +
  geom_line() +
  geom_point() +
  scale_y_continuous(limits=c(0,1), breaks=seq(0,1, by=0.1)) +
  labs(y="Probability that at least two\n people have the same birthday",
       x="Number of people") +
  theme_classic()

• Quando \(k=23\) a probabilidade de pelo menos duas pessoas com o mesmo aniversário é pouco maior do que 0,5:

birthday(23)

## [1] 0.5072972

• Quando \(k\) é maior do que 50, a probabilidade é próxima de um.

birthday(50)

## [1] 0.9703736

• Para resolver o problema dos aniversários usamos uma solução analítica para obter uma resposta exata, nem sempre é possível obter soluções analíticas para os problemas que encontramos. Quando obter uma solução analítica não for possível (ou não for nosso objetivo) podemos encontrar soluções aproximadas usando simulações pelo método de Monte Carlo.
• O método simulação de Monte Carlo, ou simplesmente simulação, se refere a uma classe de métodos estocásticos que podem ser usados para resolver problemas analíticos gerando as quantidades de interesse de forma aleatória.

• No exemplo dos aniversários podemos criar amostras com \(k\) datas de aniversário dentre as 365 possíveis e checar se todos os aniversários são em datas diferentes.
• Faremos as amostras com reposição porque para cada uma das \(k\) escolhas de data existem 365 dias igualmente prováveis independente da data escolhida anteriormente.
• Repetiremos esse exercício muitas vezes e calcularemos em quantas simulações tivemos pelo menos duas datas de aniversário iguais.

• Para fazer a simulação no R usaremos a função sample() com o argumento replace definido como TRUE.
• Além do argumento replace a função sample() tem os argumentos x, vetor do qual a amostra será retirada, size número de elementos retirados de \(x\), e prob um vetor com pesos para as probabilidades de escolher cada elemento de x.
• No nosso exemplo, como todos os elementos tem a mesma probabilidade, não usaremos o argumento prob.

sample(1:10, 3, replace = TRUE)

## [1] 5 1 6

sample(1:10, 10, replace = TRUE)

##  [1]  3  3  4  3  2  5  2 10  9  4

sample(1:10, 10, replace = FALSE)

##  [1]  4  6  1  5  7  9 10  2  3  8

• A simulação consistirá em tirar 23 números entre 1 e 365 e usar a função unique() para checar se tem números repetidos.
• Repetiremos o exercício acima mil vezes.

set.seed(4444)
k <- 23
sims <- 1000
event <- 0

for (i in 1:sims){
  days <- sample(1:365, k, replace = TRUE)
  days.unique <- unique(days)
  if (length(days.unique) < k) {
    event <- event + 1
  }
}

event/sims

## [1] 0.511

• O exercício de simulação nos deu uma probabilidade de 0,511, sabemos que a probabilidade é de 0,507. A diferença entre a probabilidade simulada e a probabilidade exata é chamada de erro de Monte Carlo.
• Embora não seja possível eliminar o erro de Monte Carlo, podemos reduzir o tamaho desse erro aumentando o número de simulações.

• Cem mil simulações.

set.seed(4444)
k <- 23
sims <- 100000
event <- 0

for (i in 1:sims){
  days <- sample(1:365, k, replace = TRUE)
  days.unique <- unique(days)
  if (length(days.unique) < k) {
    event <- event + 1
  }
}

• Cem mil simulações.

event/sims

## [1] 0.50781

• O método de simulação de Monte Carlo diz respeito a uma classe de procedimentos de repetição de procedimentos de amostragem aleatória usadas para resolver problemas analíticos de forma aproximada. Os procedimentos usados usam amostras obtidas com reposição, em que a mesma unidade pode ser escolhida repetidas vezes, e sem reposição em que cada unidade só pode ser escolhida uma vez.

• Outro método de contagem importante em probabilidade é chamado de combinações. É semelhante as permutações, porém não considera a ordem dos \(k\) elementos. Por exemplo, escolhendo dois elementos as diferentes permutações \(AB\) e \(BA\) representam a mesma combinação.
• Como a ordem dos elementos não importa, o número de combinações será sempre menor do que o número de permutações. Se tivermos de escolher 2 de 3 elementos \(A\), \(B\) e \(C\), o número de permutações será 6, \((AB, BA, AC, CA, BC, CB)\) e o número de combinações será três, \((AB, AC, BC)\).

• De fato, para calcular o número de combinações primeiro calculamos o número de permutações, \(_nP_k\) e depois dividimos por \(k!\). A divisão é feita porque para uma amostra de \(k\) elementos existem \(k!\) maneiras de arranjar os números em diferentes ordens e todos esses arranjos correspondem a uma combinação.
• No nosso exemplo, para cada dois elementos existem \(2!=2\) arranjos, por exemplo, \(AB\) e \(BA\), para \(A\) e \(B\), mas ambos contam como uma combinação.

• O número de combinações quando escolhendo \(k\) elementos de um total de \(n\) elementos é denotado por \(_nC_k\) ou \(\begin{pmatrix} n \\k \end{pmatrix}\) e calculado pela fórmula: \[_nC_k = {n \choose k} = \frac{_nP_k}{k!}=\frac{n!}{k!(n-k)!}\]

• Suponha que de um grupo de 20 pessoas, 10 homens e 10 mulheres, queiramos formar uma comissão de 5 pessoas, cada pessoa tem a mesma probabilidade de ser chamada para a comissão. Qual a probabilidade de ter pelo menos duas mulheres na comissão? \[P(\mbox{pelo menos 2 mulheres na comissão})=\\1-P(\mbox{nenhuma mulher na comissão})-\\P(\mbox{exatamente uma mulher na comissão})\]

• O número de maneiras que podemos colocar 5 de 20 pessoas em uma comissão sem considerar gênero é igual a \(_{20}C_5=15504\)
• O número de maneiras de ter uma comissão sem nenhuma mulher é igual a \(_{10}C_0 \times _{10}C_5 = 1\times 252=252\), números de maneiras de ter uma comissão sem mulheres vezes o número de maneiras de combinar cinco homens. Logo: \[P(\mbox{nenhuma mulher na comissão})=\frac{252}{15504}=0,\!016\]

• O número de maneiras de ter exatamente uma mulher na comissão é dado por \(_{10}C_1 \times _{10}C_4 = 10\times 210=2100\), Logo: \[P(\mbox{exatamente uma mulher na comissão})=\frac{2100}{15504}=0,\!135\]
• Portanto: \[P(\mbox{pelo menos 2 mulheres na comissão})=1-0,\!016-0,\!135=0,\!849\]

• Podemos simular o problema acima.

grupo <- rep(c("H", "M"), 10)

duas.mulheres <- rep(NA, 10000)

for (i in seq_along(duas.mulheres)) {
  comite <- sample(grupo, 5, replace = FALSE)
  total.mulheres = sum(comite == "M") 
  duas.mulheres[i] = total.mulheres > 1
}

sum(duas.mulheres)/length(duas.mulheres)

## [1] 0.8483

• O próximo exemplo de combinação é um tanto mais sofisiticado. Em 2009, Arnold Schwarzenegger, então governador da Califórnia, vetou um projeto de lei da assembléia.
• As primeiras letras de cada linha da mensagem formavam o que parecia ser um recado para a assembléia.

• Quando cobrado a respeito da suposta mensagem, Schwarzenegger alegou ter sido coincidência. Qual a probabilidade disso acontecer?
• O exemplo está no artigo: Philip Stark. Null and vetoed: Chance coincidence? Chance, v.23, n.4, 2009.

• Por simplicidade, suponha que o governador entregou a mensagem a um secretário que digitou o recado quebrando as linhas de forma aleatória, ou seja, as 85 palavras foram dividas em sete quebras de linha aleatórias, cada uma com pelo menos uma palavra.
• Também suponha que as palavras não podem ser quebradas, as maneiras de quebrar as linhas têm a mesma probabilidade e o número total de linhas é fixo e igual a sete.

• A primeira coisa a fazer é calcular como 85 palavras podem ser divididas em sete linhas, ou seja, como inserir seis quebras de linha. Uma quebra de linha pode ser inserida depois da primeira palavra, depois da segunda, da terceira… ou antes da 85o. Logo, existem 84 maneiras de inserir 6 quebras de linha.

• Como a ordem das quebras de linha não é relevante (a ordem das palavras é relevante!) usaremos combinações e não permutas. Para calcular o número de combinações usaremos a função choose() do R, a função lchoose() faz o cálculo com logaritmos.

choose(84,6)

## [1] 406481544

• Existem mais de 400 milhões de formas de colocar as quebras de linhas!

• Quantas são as formas de quebrar as linha de modo a aparecer a suposta mensagem?
• Para a segunda linha começar com u a quebra tinha de estar em unnecessary, só uma possibilidade. Para a terceira linha começar com c a quebra tinha de estar em come, consideration ou care, três possibilidades. Para o k na quarta linha a quebra tinha de estar em kicks, uma possibilidade. Para a quebra do y são duas possibilidades: Yet ou year. Para o o também duas possibilidades: overwhelmingly ou of. Finalmente, para o u na última linha há uma possibilidade: unnecessary*. Logo, existem \(1\times 3 \times 1 \times 2 \times 2 \times 1 = 12\) maneiras da suposta mensagem aparecer.

• A probabilidade de ter sido obra do acaso é \(\frac{12}{_{84}C_6}=\frac{12}{406.481.544}\) que é aproximadamente de uma em 34 milhões!

• Probabilidade condicional está relacionada a mudança na probabilidade de um evento depois que observamos outros eventos.
• A probabilidade condicional permite que evidências observadas sejam levadas em consideração quando do cálculo da probabilidade de um evento ocorrer.

• A robabilidade condicional de um evento \(A\) ocorrer dado que o evento \(B\) ocorreu é denotada por \(P(A|B)\) e definida por: \[P(A|B)=\frac{P(A\mbox{ e }B)}{P(B)}\]
• Na equação acima, \(P(A\mbox{ e }B)\) é a probabilidade conjunta dos dois eventos ocorrerem e \(P(B)\) é a probabilidade marginal do evento \(B\).

• Arrumando os termos da equação podemos obter a seguinte regra de multiplicação: \[P(A\mbox{ e }B) = P(A|B)\times P(B)=P(B|A)\times P(A)\]
• Com esta regra podemos encontrar uma forma alternativa para lei da probabilidade total: \[P(A)=P(A|B)\times P(B) + P(A|B^c)\times P(B^c)\]

• Para ajudar a entender a importância da probabilidade condicional, considere um exemplo onde dois casais estão esperando gêmeos. Um casal fez um exame de ultrassom e conseguir determinar que um dos bebês é menino. O outro casal não sabia o sexo dos bebês até o nascimento, quando viu que o primeiro a nascer foi um menino.
• Qual a probabilidade dos dois bebês serem meninos? Essa probabilidade é a mesma para os dois casais? A diferença entre os casais é que um casal sabe que pelo menos um dos bebês é menino e o outro sabe que o primeiro bebê é menino.

• Comecemos registrando os quatro resultados no espaço amostral. Denote o sexo do bebê por \(G\) (girl) ou \(B\) (boy), o espaço amostral é dado por: \(\Omega = \{GG, GB, BG, BB \}\).

• Para o primeiro casal a probabilidade de interesse é: \[P(BB|\mbox{pelo menos um é menino})=\frac{P(BB\mbox{ e pelo menos um é menino})}{P(\mbox{pelo menos um é menino})}\\=\frac{P(BB\mbox{ e }\{BB\mbox{ ou }BG\mbox{ ou }GB\})}{P(BB\mbox{ ou }BG\mbox{ ou }GB)}\\= \frac{P(BB)}{P(BB\mbox{ ou }BG\mbox{ ou }GB)}=\frac{1/4}{3/4}=\frac{1}{3}\]

• Para o segundo casal a probabilidade de interesse é: \[P(BB|\mbox{o mais velho é menino})=\frac{P(BB\mbox{ e mais velho é menino})}{P(\mbox{mais velho é menino})}\\=\frac{P(BB\mbox{ e }\{BB\mbox{ ou }BG\})}{P(BB\mbox{ ou }BG)}\\=\frac{P(BB)}{P(BB\mbox{ ou }BG)}=\frac{1/4}{1/2}=\frac{1}{2}\]

• O exemplo nos mostra a relevância da informação que condiciona a probabilidade. Em um caso a probabilidade é condicionada a saber que pelo menos um é menino, no outro caso é condicionada a saber que o primeiro é menino.
• No primeiro caso temos três possibilidades, \(\{GB, BG, BB \}\), no segundo caso, sabendo que o primeiro bebê é menino, temos duas possibilidades, \(\{BG, BB \}\).

• Probabilidades condicionais podem ser usadas para descrever caracterísiticas de uma população, por exmeplo, se 10% dos eleitores são afro-americanos, então podemos escrever \(P(\mbox{black})=0,\!1\). Podemos interpretar esse número dizendo que se escolhermos aleatoriamente um eleitor a probabilidade dele ser afro-americano é de 10%.
• De modo semelhante, \(P(\mbox{black}|\mbox{hispanic ou black})\) representa a proporção de afro-americanos entre eleitores afro-americanos e hispânicos.

• O arquivo FLVoters tem uma amostra aleatória de dez mil eleitores da Flórida.
• Usaremos esse arquivo para ilustrar o uso de probabilidade condicional.

• Variáveis:
  • surname: sobrenome
  • county: id do condado onde o eleitor mora
  • VTD: id do distrito onde o eleitor mora
  • age: idade
  • gender: gênero
  • race: raça

FLVoters <- read.csv("FLVoters.csv")
glimpse(FLVoters)

## Rows: 10,000
## Columns: 6
## $ surname <chr> "PIEDRA", "LYNCH", "CHESTER", "LATHROP", "HUMMEL", "CHRISTISON…
## $ county  <int> 115, 115, 115, 115, 115, 115, 115, 115, 1, 1, 115, 115, 115, 1…
## $ VTD     <int> 66, 13, 103, 80, 8, 55, 84, 48, 41, 39, 26, 45, 11, 48, 22, 88…
## $ age     <int> 58, 51, 63, 54, 77, 49, 77, 34, 56, 60, 44, 45, 80, 83, 88, 55…
## $ gender  <chr> "f", "m", "m", "m", "f", "m", "f", "f", "f", "m", "m", "f", "m…
## $ race    <chr> "white", "white", NA, "white", "white", "white", "white", "whi…

• Para começar vamos remover as observações com informações ausentes.

FLVoters <- FLVoters %>%
  na.omit()

dim(FLVoters)

## [1] 9113    6

• Perdemos 887 observações.

• Para o exemplo, vamos tratar essa amostra de 9.113 eleitores como se fosse a população de interesse. Vamos calcular a probabilidade marginal de cada raça. Para isso vamos dividir o número de observações que declararam a raça pelo total de obsevações.

margin_race <- FLVoters %>%
  count(race) %>%
  mutate(prop = n/sum(n))

margin_race

##       race    n        prop
## 1    asian  175 0.019203336
## 2    black 1194 0.131021617
## 3 hispanic 1192 0.130802151
## 4   native   29 0.003182267
## 5    other  310 0.034017338
## 6    white 6213 0.681773291

• Distribuição marginal por gênero:

margin_gender <- FLVoters %>%
  count(gender) %>%
  mutate(prop = n/sum(n))

margin_gender

##   gender    n      prop
## 1      f 4883 0.5358279
## 2      m 4230 0.4641721

• Distribuição condicional de raça dado que o gênero é feminino.

margin_race_f <- FLVoters %>%
  filter(gender == "f") %>%
  count(race) %>%
  mutate(prop = n/sum(n))

margin_race_f

##       race    n        prop
## 1    asian   83 0.016997747
## 2    black  678 0.138849068
## 3 hispanic  666 0.136391563
## 4   native   17 0.003481466
## 5    other  158 0.032357157
## 6    white 3281 0.671922998

• Distribuição condicional de raça dado que o gênero é masculino.

margin_race_m <- FLVoters %>%
  filter(gender == "m") %>%
  count(race) %>%
  mutate(prop = n/sum(n))

margin_race_m

##       race    n        prop
## 1    asian   92 0.021749409
## 2    black  516 0.121985816
## 3 hispanic  526 0.124349882
## 4   native   12 0.002836879
## 5    other  152 0.035933806
## 6    white 2932 0.693144208

• Por fim, podemos calcular a distribuição conjunta de raça e gênero.

joint_p <- FLVoters %>%
  count(gender, race) %>%
  mutate(prop = n/sum(n))

##    gender     race    n        prop
## 1       f    asian   83 0.009107868
## 2       f    black  678 0.074399210
## 3       f hispanic  666 0.073082410
## 4       f   native   17 0.001865467
## 5       f    other  158 0.017337869
## 6       f    white 3281 0.360035115
## 7       m    asian   92 0.010095468
## 8       m    black  516 0.056622408
## 9       m hispanic  526 0.057719741
## 10      m   native   12 0.001316800
## 11      m    other  152 0.016679469
## 12      m    white 2932 0.321738176

• Da distribuição conjunta podemos calcular a distribuição marginal e a condicional.
• Por exemplo, a probabilidade de ser afro-americano é dada por: \[P(\mbox{black}) = P(\mbox{black e female}) + P(\mbox{black e male})\]
• De fato, 0,074 + 0,057 = 0,131 que é o que diz a distribuição marginal.

• Podemos usar o R para obter a distribuição marginal a partir da conjunta. Para isso usaremos a função pivot_wider() para criar uma coluna para cada distribuição conjunta e depois somamos as duas.

joint_p %>%
  select(-n) %>%
  pivot_wider(names_from = gender, values_from = prop)

## # A tibble: 6 × 3
##   race           f       m
##   <chr>      <dbl>   <dbl>
## 1 asian    0.00911 0.0101 
## 2 black    0.0744  0.0566 
## 3 hispanic 0.0731  0.0577 
## 4 native   0.00187 0.00132
## 5 other    0.0173  0.0167 
## 6 white    0.360   0.322

joint_p %>%
  select(-n) %>%
  pivot_wider(names_from = gender, values_from = prop) %>%
  mutate(total_prop = f + m)

## # A tibble: 6 × 4
##   race           f       m total_prop
##   <chr>      <dbl>   <dbl>      <dbl>
## 1 asian    0.00911 0.0101     0.0192 
## 2 black    0.0744  0.0566     0.131  
## 3 hispanic 0.0731  0.0577     0.131  
## 4 native   0.00187 0.00132    0.00318
## 5 other    0.0173  0.0167     0.0340 
## 6 white    0.360   0.322      0.682

• De forma semelhante, podemos obter a distribuição marginal de gênero, mas como existem seis raças, a soma terá seis termos.
• As funções rowSums(), que soma as linhas de um data.frame, across(), que permite selecionar colunas com a mesma semântica da função select(), e where(), que seleciona as variáveis para as quais a condição estabelecida é verdadeira nos ajudarão a fazer a soma.

joint_p %>%
  select(-n) %>%
  pivot_wider(names_from = race, values_from = prop) %>%
  mutate(total_prop = rowSums(across(where(is.numeric))))

## # A tibble: 2 × 8
##   gender   asian  black hispanic  native  other white total_prop
##   <chr>    <dbl>  <dbl>    <dbl>   <dbl>  <dbl> <dbl>      <dbl>
## 1 f      0.00911 0.0744   0.0731 0.00187 0.0173 0.360      0.536
## 2 m      0.0101  0.0566   0.0577 0.00132 0.0167 0.322      0.464

• A dsitribuição condicional é obtida dividindo a distribuição conjunta pela marginal. \[P(\mbox{black}|\mbox{female}) = \frac{P(\mbox{black e female})}{P(\mbox{female})}\\=\frac{0,\!074}{0,\!536}=0,\!139\]

• O conceito de probabilidade condicional pode ser usado em conjuntos com mais de dois eventos. Para os eventos \(A\), \(B\) e \(C\), com probabilidade conjunta dada por \(P(\mbox{A e B e C})\) temos dois tipo de probabilidade condicional: a probabilidade de dois eventos condicionada ao terceiro, \(P(\mbox{A e B}|C)\), e a probabilidade de um evento condicionada aos outros dois, \(P(A|\mbox{B e C})\). \[P(\mbox{A e B}|C) = \frac{P(\mbox{A e B e C})}{P(C)}\] \[P(A|\mbox{B e C})=\frac{P(\mbox{A e B e C})}{P(\mbox{B e C})}=\frac{P(\mbox{A e B|C})\times P(C)}{P(\mbox{B e C})}\\=\frac{P(\mbox{A e B|C})\times P(C)}{P(\mbox{B|C})\times P(C)}=\frac{P(\mbox{A e B|C})}{P(\mbox{B|C})}\]

• Para ilustrar as probabilidade condicionais acima vamos criar uma variável com grupos de idade. Podemos fazer isso com a função case_when(), como fizemos em outros exemplos, ou com a função cut().
• O argumento right da função cut() diz para o R que o valor da fronteira superior deve ficar no grupo inferior, ou seja, o intervalo é fechado a direita.

FLVoters <- FLVoters %>%
  mutate(age_group = cut(age, breaks = c(0,20,40,60,Inf),
                         right = TRUE,
                         labels = c("<= 20", "20-40", "40-60", "> 60")))

• A função cut() cria a variável como fator.

glimpse(FLVoters)

## Rows: 9,113
## Columns: 7
## $ surname   <chr> "PIEDRA", "LYNCH", "LATHROP", "HUMMEL", "CHRISTISON", "HOMAN…
## $ county    <int> 115, 115, 115, 115, 115, 115, 115, 1, 1, 115, 115, 115, 115,…
## $ VTD       <int> 66, 13, 80, 8, 55, 84, 48, 41, 39, 26, 45, 11, 48, 88, 25, 8…
## $ age       <int> 58, 51, 54, 77, 49, 77, 34, 56, 60, 44, 45, 80, 83, 55, 33, …
## $ gender    <chr> "f", "m", "m", "f", "m", "f", "f", "f", "m", "m", "f", "m", …
## $ race      <chr> "white", "white", "white", "white", "white", "white", "white…
## $ age_group <fct> 40-60, 40-60, 40-60, > 60, 40-60, > 60, 20-40, 40-60, 40-60,…

• Distribuição conjunta de raça, gênero e idade

joint3 <- FLVoters %>%
  count(race, age_group, gender) %>%
  mutate(prop = n/sum(n))

-Qual a probabilidade de ser mulher e afro-americana dado que tem mais de 60 anos? - Para responder essa pergunta precisamos dividir a probabilidade conjunta, que está em joint3, pela probabilidade marginal de ter mais de 60.

• Probabilidade marginal por grupos de idade

margin_age <- FLVoters %>%
  count(age_group) %>%
  mutate(margin_age = n/sum(n)) %>%
  select(-n)

margin_age

##   age_group margin_age
## 1     <= 20 0.01766707
## 2     20-40 0.27093164
## 3     40-60 0.36047405
## 4      > 60 0.35092725

joint3 <- left_join(joint3, margin_age, by="age_group") %>%
  mutate(con_prob_age = prop/margin_age)

filter(joint3, race == "black", gender == "f", age_group == "> 60") %>%
  select(con_prob_age)

##   con_prob_age
## 1   0.05378361

• Para a probabilidade de ser afro-americano dado que é mulher e tem mais de 60 anos devemos dividir a distribuição conjunta das três variáveis pela distribuição conjunta de idade e gênero. \[\frac{P(\mbox{black e 60+ e female})}{P(\mbox{60+ e female})}\]

joint2 <- FLVoters %>%
  count(age_group, gender) %>%
  mutate(prob_age_gender = n/sum(n)) %>%
  select(-n)

joint2

##   age_group gender prob_age_gender
## 1     <= 20      f     0.009656535
## 2     <= 20      m     0.008010534
## 3     20-40      f     0.143092286
## 4     20-40      m     0.127839350
## 5     40-60      f     0.189838692
## 6     40-60      m     0.170635356
## 7      > 60      f     0.193240426
## 8      > 60      m     0.157686821

joint3 <- left_join(joint3, joint2, by=c("age_group", "gender")) %>%
  mutate(con_prob_race = prop/prob_age_gender)

filter(joint3, gender == "f", age_group == "> 60", race == "black") %>%
  select(con_prob_race)

##   con_prob_race
## 1    0.09767178

• Um importante conceito ligado à probabilidade condicional é o de independência. De forma intuitiva, independência entre dois eventos significa que informações a respeito da ocorrência do primeiro evento não nos oferece informações adicionais a respeito da ocorrência do segundo evento.

• Se os eventos \(A\) e \(B\) são independentes um do outro, a probabilidade condicional de \(A\) dado \(B\) não difere da probabilidade marginal de \(A\), vale o mesmo para \(B\) condicional a \(A\). \[P(A|B)=P(A) \qquad \mbox{e} \qquad P(B|A)=P(B)\]

• Formalmente podemos definir independência a partir da regra de multiplicação que vimos anteriormente segundo a qual: \[P(A\mbox{ e }B) = P(A|B)\times P(B)=P(B|A)\times P(A)\]

• Os eventos \(A\) e \(B\) são independentes se e somente a distribuição de probabilidade conjunta é igual ao produto das distribuições marginais: \[P(A\mbox{ e }B) = P(A)\times P(B)\]

• Como exemplo, vamos avaliar se raça e gênero são independentes na amostra de eleitores da Flórida que estamos utilizando. Naturalmente não esperamos encontrar uma relação exatamente independente no sentido da definição formal, vamos avaliar se, por exemplo, a proporção de mulheres é maior do que a esperada em determinados grupos raciais.

• Se valer independência devemos ter que \(P(\mbox{black e female}) = P(\mbox{black})\times P(\mbox{female})\), \(P(\mbox{white e male}) = P(\mbox{white})\times P(\mbox{male})\), e assim por diante.
• Vamos usar um gráfico de dispersão para comparar a distribuição conjunta com o produto das marginais. Quanto mais perto da linha de 45o os pontos estiverem, mais próximo de independentes serão os eventos.

margin_race <- margin_race %>%
  select(race, prob_race = prop)

margin_gender <- margin_gender %>%
  select(gender, prob_gender = prop)

race_gender_indep <- tidyr::crossing(margin_race, margin_gender) %>%
  left_join(select(joint_p, gender, race, prob_joint=prop),
            by = c("gender", "race")) %>%
  mutate(prob_indp = prob_race * prob_gender)

race_gender_indep %>%
  filter(gender == "f") %>%
  ggplot(aes(prob_indp, prob_joint)) +
  geom_abline(intercept=0, slope=1, color = "black", linewidth = 0.5) +
  geom_point(shape=1, size=2) +
  coord_fixed() +
  labs(x= "P(race)*P(gender)",
       y="P(race e gender)") +
  theme_classic()

• O conceito de independência pode ser aplicado a mais de dois eventos. Por exemplo, se temos os eventos \(A\), \(B\) e \(C\) a independência conjunta desses eventos implicas que a probabilidade conjunta pode ser escrita como o produtos das probabilidades marginais.

\[P(\mbox{A e B e C})=P(A)\times P(B)\times P(C)\] - A independência pode ser definida em termos condicionais. \[P(\mbox{A e B|C})=P(\mbox{A|C})\times P(\mbox{B|C})\]

• Se \(A\), \(B\) e \(C\) são independentes, então \(A\) e \(B\) são independentes. Para ver isso, note que: \[P(A\mbox{ e }B) = P(A\mbox{ e }B\mbox{ e } C) + P(A\mbox{ e }B\mbox{ e }C^c)\\=P(A)P(B)P(C)+P(A)P(B)P(C^c)\\=P(A)P(B)[P(C)+P(C^c)]=P(A)P(B)\]

• Independência conjunta implica em independência condicional, mas o contrário não é verdade. \[P(A\mbox{ e }B|C)=\frac{P(A\mbox{ e }B\mbox{ e }C)}{P(C)}=\frac{P(A)P(B)P(C)}{P(C)}=P(A|C)P(B|C)\]

• Para ilustrar esses conceitos vamos comparar as distribuições conjuntas em joint3 com os produtos das probabilidades marginais em margin_race, margin_age e margin_gender.
• Depois vamos examinar a dependência condicional entre raça e idade.

margin_age <- margin_age %>% 
  select(age_group, prob_age=margin_age)

joint_indep <- tidyr::crossing(margin_age, margin_gender, margin_race) %>%
  mutate(indep_prob = prob_race * prob_age * prob_gender) %>%
  left_join(select(joint3, race, age_group, gender, prob_joint=prop),
          by=c("gender", "age_group", "race"))

• Gráfico para comparação

joint_indep %>%
  filter(age_group == "> 60", gender == "f") %>%
  ggplot(aes(prob_joint, indep_prob)) +
  geom_abline(intercept=0, slope = 1, color="black", linewidth=0.5) +
  geom_point(shape=1, size=2) +
  coord_fixed() +
  labs(title = "Joint independence",
       x="P(race e 60+ e female)",
       y="P(race)*P(60+)*P(female)") +
  theme_classic()

#raça e idade dado gênero
cond_prob_gender <- left_join(select(joint3, race, age_group, gender, joint_prob = prop),
                              margin_gender, by=c("gender")) %>%
  mutate(cond_prob_gender = joint_prob/prob_gender)

#raça dado gênero
cond_prob_race_gender <- left_join(select(joint_p, race, gender, prob_race_gender=prop),
                                   margin_gender, by=c("gender")) %>%
  mutate(cond_prob_race_gender = prob_race_gender/prob_gender) %>%
  select(race, gender, cond_prob_race_gender)

#idade dado gênero
cond_prob_age_gender <- left_join(select(joint2, age_group, gender, prob_age_gender), 
                                  margin_gender, by=c("gender")) %>%
  mutate(cond_prob_age_gender = prob_age_gender/prob_gender) %>%
  select(age_group, gender, cond_prob_age_gender)

#probabilidade conjunta de raça e idade
indep_cond_gender <- full_join(cond_prob_race_gender, cond_prob_age_gender,
                               by = c("gender")) %>%
  mutate(indep_prob = cond_prob_race_gender * cond_prob_age_gender)

#juntar probabilidades conjuntas com as condicionais
master <- left_join(cond_prob_gender, indep_cond_gender,
                    by = c("gender", "age_group", "race"))

• Gráfico

master %>%
  filter(age_group == "> 60", gender == "f") %>%
  ggplot(aes(cond_prob_gender, indep_prob)) +
  geom_abline(intercept = 0, slope = 1, color="black", linewidth=0.5) +
  geom_point(shape=1, size=2) +
  coord_fixed() +
  labs(title = "Marginal independence",
       x = "P(race and 60+ | female)",
       y="P(race | female) * P(60+ | female)") +
  theme_classic()

• Você quer trocar de porta?
• Nos anos 80 Sérgio Mallandro apresentava um programa que tinha um quadro chamado de Porta dos Desesperados. O jogo consistia de três portas, uma tinha um prêmio e as outras tinham pegadinhas (por exemplo, um sujeito fantasiado de gorila).
• O participante escolhia uma porta. Sérgio Mallandro abria outra porta que não tinha o prêmio e perguntava se o paticipante queria trocar de porta.
• E aí? Troca ou não troca?

• De saída pode parecer que tanto faz trocar ou não de porta. Afinal, depois do Mallandro abrir uma porta, sobram duas e qualquer delas pode ter o prêmio. Logo, trocando ou não, a chance de ganhar seria de 50%, mas…
• Vamos pensar com cuidado e usar o que aprendemos nesta unidade.

• Considere a estratégia de não trocar, nesse caso o resultado está determinado sem importar o Mallandro faça. Logo, a probabilidade de ganhar o prêmio é de 1/3 e a de receber um pegadinha é de 2/3.

• Agora considere a estratégia de trocar. Existem dois cenários:
  • No primeiro você acertou a porta na primeira escolha, nesse caso trocar é um mau negócio. A probabilidade desse caso é de 1/3.
  • No segundo você errou a porta na primeira escolha, nesse caso trocar é um bom negócio. A probabilidade de errar a porta na primeira escolha é de 2/3. Como o Mallandro abiu uma porta sem prêmio e a escolhida por você também não tem o prêmio se você trocar com certeza ganhará o prêmio.

• Em resumo. Com a estratégia não trocar você ganha se a escolha inicial for correta, probabilidade de 1/3, com a estratégia de trocar você ganha se a escolha inicial for a porta errada, probabilidade de 2/3. A probabilidade de ganhar com a estratégia trocar é duas vezes maior do que com a estratégia não trocar!

• Estranho? Desconfiado? Façamos a conta!
• Proabilidade de ganhar o prêmio: \[P(\mbox{prêmio})=P(\mbox{prêmio}|\mbox{prêmio na primeira})P(\mbox{prêmio na primeira})\\+P(\mbox{prêmio}|\mbox{pegadinha na primeira})P(\mbox{pegadinha na primeira})\]

• Seja qual for a estratégia, vale que:

\[P(\mbox{prêmio na primeira}) = \frac{1}{3}\] \[P(\mbox{pegadinha na primeira}) = \frac{2}{3}\]

• Se a estratégia for não trocar, temos: \[P(\mbox{prêmio}|\mbox{prêmio na primeira})=1\] \[P(\mbox{prêmio}|\mbox{pegadinha na primeira})=0\]

• Logo, com a estratégia de não trocar: \[P(\mbox{prêmio}) = 1\times \frac{1}{3} + 0\times \frac{2}{3} = \frac{1}{3}\]

• Se a estratégia for trocar, temos: \[P(\mbox{prêmio}|\mbox{prêmio na primeira})=0\] \[P(\mbox{prêmio}|\mbox{pegadinha na primeira})=1\]

• Logo, com a estratégia de não trocar: \[P(\mbox{prêmio}) = 0\times \frac{1}{3} + 1\times \frac{2}{3} = \frac{2}{3}\]

• Ainda descrente? Vamos fazer uma simulação de Monte Carlo!

sims <- 1000
portas <- c("pegadinha", "pegadinha", "premio")

result.trocar <- result.naotrocar <- rep(NA, sims)

for (i in 1:sims){
  #escolha inicial
  primeira.escolha <- sample(1:3, size=1)
  result.naotrocar[i] <- portas[primeira.escolha]
  restante <- portas[-primeira.escolha]
  #Mallandro abre uma porta com pegadinha
  if(portas[primeira.escolha]=="premio") #escolhe uma das duas com pegadinha
    mallandro <- sample(1:2, size=1)
  else #só sobrou uma com pegadinha
    mallandro <- (1:2)[restante == "pegadinha"]
  result.trocar[i] <- restante[-mallandro]
}

mean(result.naotrocar == "premio")

## [1] 0.308

mean(result.trocar == "premio")

## [1] 0.692

• Conforme sabíamos a estratégia não trocar deu certo cerca de 1/3 das vezes e a estratégia trocar deu certo cerca de 2/3 das vezes.
• No livro o exemplo aparece como problema Monty Hall em referência ao apresentados do programa nos Estados Unidos.

• No começo desta unidade falamos de diferentes interpretações para probabilidade: frequentista e bayesiana.
• Na perspectiva bayesiana, a subjetiva, é natural perguntar como fazemos para atualizar nossas crenças iniciais após observar os dados.

• A regra de Bayes é dada pela equação: \[P(A|B)=\frac{P(B|A)P(A)}{P(B)}=\frac{P(B|A)P(A)}{P(B|A)P(A)+P(B|A^c)P(A^c)}\]
• Na equação, \(P(A)\) é chamada probabilidade anterior ou a priori (prior) e reflete a crença inicial sobre a probabilidade do evento \(A\) ocorrer. Depois de observar os dados, representados por \(B\), a crença é atualizada para \(P(A|B)\), que é chamada de probabilidade posterior (posterior).

• A regra de Bayes é válida independente da forma como a probabilidade é interpretada, isso é verdade porque a regra decorre da definição de probabilidade condicional e da regra da probabilidade total. \[P(A|B)=\frac{P(A\mbox{ e }B)}{P(B)}=\frac{P(B|A)P(A)}{P(B)}\]

• Uma aplicação conhecida da regra de Bayes é a interpretação de testes de diagnóstico em medicina, esses testes podem dar falsos negativos e falsos positivos.
• Considere o seguinte exemplo. Uma mulher grávida de 35 anos fica sabendo que 1 em 378 mulheres da idade dela tem filhos com síndrome de Down (DS).

• Um ultrassom no primeiro trimestre de gravidez a coloca no grupo de risco. De cada 100 casos de DS, em 86 a mãe recebe o resultado de alto risco, nos outros 14 não recebe resultado de alto risco. A probabilidade de uma mãe que não é de alto risco receber um diagnóstico de alto risco é de 1 em 20.
• Dado que mãe recebeu o diagnóstico de alto risco, qual a probabilidade do filho ter DS? Qual seria a probabilidade se o resultado fosse negativo (não deu alto risco)?

• Para resolver esse problema começamos definindo a prior, sem teste algum a probabilidade do filho ter DS é de 1/378, aproximadamente 0,003.
• O exame de ultrassom retorna um resultado de alto risco 86% das vezes que o filho tem DS, \(P(HR|DS)=0,\!86\). Por outro lado, em 5% dos casos (1/20) o teste da alto risco para quem não é de alto risco, \(P(HR|\mbox{not }DS)=0,\!05\).

• Com essas informações e a regra de Bayes podemos calcular a probabilidade do filho ter DS dado que o teste deu alto risco: \[P(DS|HR)=\frac{P(HR|DS)P(DS)}{P(HR|DS)P(DS)+P(HR|\mbox{not }DS)P(\mbox{not }DS)}\\=\frac{0,\!86 \times 1/378}{0,\!86\times 1/378+0,\!05\times 377/378}=0,\!04\]

• De forma semelahnte podemos calcular a probabilidade do filho ter DS se a mãe recebeu um negativo no ultrassom. \[P(DS|\mbox{not }HR)=\\ \frac{P(\mbox{not }HR|DS)P(DS)}{P(\mbox{not }HR|DS)P(DS)+P(\mbox{not }HR|\mbox{not }DS)P(\mbox{not }DS)}\\=\frac{0,\!14 \times 1/378}{0,\!14\times 1/378+0,\!95\times 377/378}=0,\!0004\]

• Mesmo quando a mãe recebe um diagnóstico de alto risco, a probabilidade do filho ter DS é baixa. Isso ocorre porque DS é uma condição rara, como mostra o pequeno valor para prior.

• Para um outro exemplo de regra de Bayes considere uma urna com duas bolas. Sabemos que uma das bolas é preta, não sabemos se a outra é branca ou preta. Não existem bolas de outras cores.
• Logo, existem apenas duas hipóteses: PP (preta e preta), PB (preta e branca). Antes de obervar a cor de uma das bolas, você acredita que a probabilidade de ser PP é igual a probabilidade de ser BB, ou seja, \(P(PP)=P(PB)=0,\!5\).
• Uma bola é retirada da urna e observamos que a cor é preta.

• Dado que a primeira bola é preta, vamos chamar esse evento de \(1P\) qual a probabilidade da segunda ser preta? \[P(PP|1P)=\frac{P(1P|PP)P(PP)}{P(1P)}\\=\frac{P(1P|PP)P(PP)}{P(1P|PP)P(PP)+P(1P|PB)P(PB)}\\=\frac{1\times0,\!5}{1\times0,\!5+0,\!5\times 0,\!5}=\frac{0,\!5}{0,\!75}=\frac{2}{3}\]

• \(P(1P|PP)\): probabilidade da primeira ser preta dado que as duas são pretas, igual a 1
• \(P(PP)\): prior para PP, igual a 0,5
• \(P(1P|PB)\): probabilidade da primeira ser preta dado que uma é branca e outra preta, igual a 0,5
• \(P(PB)\): prior para PB, igual a 0,5

• Qual a probabilidade da segunda ser branca? \[P(PB|1P)=\frac{P(1P|PB)P(PB)}{P(1P)}\\=\frac{P(1P|PB)P(PB)}{P(1P|PP)P(PP)+P(1P|PB)P(PB)}\\=\frac{0,\!5\times0,\!5}{1\times0,\!5+0,\!5\times 0,\!5}=\frac{0,\!25}{0,\!75}=\frac{1}{3}\]

• \(P(1P|PB)\): probabilidade da primeira ser preta dado que as duas são pretas, igual a 0,5
• \(P(PP)\): prior para PP, igual a 0,5
• \(P(1P|PP)\): probabilidade da primeira ser preta dado que uma é branca e outra preta, igual a 0,5
• \(P(PB)\): prior para PB, igual a 0,5

• Agora podemos partir para exemplos mais sofisiticados do uso da regra de Bayes.
• O US Census Bureau divulga uma lista com sobrenomes comuns e as frequências desses sobrenomes, os dados incluem mais de 150.000 sobrenomes que apareceram pelo menos 100 vezes.

• Também está disponível a frequência relativa de cada raça em cada sobrenome, são seis categorias de raça: non-Hispanic White, non-Hispanic African American, non-Hispanic Asian and Pacific Islander, Hispanic origin, non-Hispanic American Indian and Alaska Native, e non-Hispanic of two or more races. Para o exemplo juntaremos as duas últimas categorias em non-Hispanic others.

• O exercício consiste inicialmente em prever a raça usando informações a respeito do sobrenome, depois usaremos outra base de dados para prever a raça usando também o endereço.
• O exemplo segue: Kosuke Imai e Kabir Khanna. Improving ecological inference by predicting individual ethnicity from voter registration record. Political Analysisi, v.24, n.2, 2016.

• Para acompanhar o exemplo é bom ter em mente que \(P(\mbox{race})\), \(P(\mbox{race}|\mbox{surname})\) e \(P(\mbox{race e surname})\) são exemplos de representar as distribuições marginal, condicional e conjunta.
• Os dados estão em cnames.

• Variáveis:
  • surname: sobrenome
  • count: número de indivíduos com o sobrenome específico
  • pctwhite: percentagem de brancos não-hispânicos entre os que possuem o sobrenome específico
  • pctblack: percentagem de afro-americanos não-hispânicos entre os que possuem o sobrenome específico
  • pctapi: percentagem de asiáticos e das ilhas do pacífico não-hispânicos entre os que possuem o sobrenome específico
  • pcthispanic: percentagem de hispânicos entre os que possuem o sobrenome específico
  • pctothers: percentagem de outros entre os que possuem o sobrenome específico

cnames <- read.csv("cnames.csv")
glimpse(cnames)

## Rows: 151,671
## Columns: 7
## $ surname     <chr> "SMITH", "JOHNSON", "WILLIAMS", "BROWN", "JONES", "MILLER"…
## $ count       <int> 2376206, 1857160, 1534042, 1380145, 1362755, 1127803, 1072…
## $ pctwhite    <dbl> 73.342666, 61.550000, 48.520000, 60.716072, 57.690000, 85.…
## $ pctblack    <dbl> 22.217778, 33.800000, 46.720000, 34.543454, 37.730000, 10.…
## $ pctapi      <dbl> 0.399960, 0.420000, 0.370000, 0.410041, 0.350000, 0.419958…
## $ pcthispanic <dbl> 1.559844, 1.500000, 1.600000, 1.640164, 1.440000, 1.429857…
## $ pctothers   <dbl> 2.479752, 2.730000, 2.790000, 2.690269, 2.790000, 1.939806…

• Na base de dados constam 151.671 sobrenomes. Para cada sobrenome, temos a probbailidade de pertencer a um determinado grupo racial, \(P(\mbox{race}|\mbox{surname})\).
• Para manter anonimato o US Census Bureau não disponibiliza percentagens muito pequenas para um dado sobrenome, isso deixa os dados com alguns NAs. Para substituir os NAs os autores subtrairam a soma das percentagens disponíveis de 100% e dividiram o restante igualmente entre os grupos raciais com dados ausentes.

• O primeiro exercício consiste em usar as probabilidades condicionais para classficar os eleitores por raça. A validação das previsões será feita usando os dados de eleitores da Flórida que estamos usando de outros exemplos.

• Inicialmente vamos identificar a raça mais provável para cada sobrenome.

most_likely_race <- cnames %>%
  select(-count) %>%
  pivot_longer(names_to = "race_pred", values_to = "highest_pct", -surname) %>%
  mutate(race_pred = str_replace(race_pred, "pct", "")) %>%
  group_by(surname) %>%
  filter(row_number(desc(highest_pct))==1) %>%
  ungroup()

cnames <- full_join(cnames, most_likely_race, by=c("surname"))

• O próximo passo é juntar o data.frame cnames com o FLVoters, para isso usaremos a função inner_join() que só aproveita as entradas comuns aos dois data.frames.

FLVotersJoin <- FLVoters %>%
  inner_join(cnames, by="surname")

• Checar se a classificação de raças é a mesma nas colunas para raça e raça prevista.

unique(FLVotersJoin$race)

## [1] "white"    "other"    "black"    "asian"    "hispanic" "native"

unique(FLVotersJoin$race_pred)

## [1] "hispanic" "white"    "black"    "others"   "api"

• Reclassficiar para que fiquem consistentes

FLVotersJoin <- FLVotersJoin %>%
  mutate(race = recode(race, "native"="other"),
         race_pred = recode(race_pred, "api"="asian", "others"="other"))

unique(FLVotersJoin$race)

## [1] "white"    "other"    "black"    "asian"    "hispanic"

unique(FLVotersJoin$race_pred)

## [1] "hispanic" "white"    "black"    "other"    "asian"

• Avaliar a taxa de acerto (true positive rate)

race_tp <- FLVotersJoin %>%
  group_by(race) %>%
  summarize(tp = mean(race == race_pred)) %>%
  arrange(desc(tp))

race_tp

## # A tibble: 5 × 2
##   race          tp
##   <chr>      <dbl>
## 1 white    0.950  
## 2 hispanic 0.847  
## 3 asian    0.564  
## 4 black    0.160  
## 5 other    0.00361

• Agora os falsos positivos

race_fp <- FLVotersJoin %>%
  group_by(race) %>%
  summarize(fp = mean(race != race_pred)) %>%
  arrange(desc(fp))

race_fp

## # A tibble: 5 × 2
##   race         fp
##   <chr>     <dbl>
## 1 other    0.996 
## 2 black    0.840 
## 3 asian    0.436 
## 4 hispanic 0.153 
## 5 white    0.0498

• O modelo apenas com o sobrenome conseguiu classificar corretamente 84,7% dos eleitores hispânicos e 95% dos eleitores brancos. Nas outras raças o modelo não se saiu bem, por exemplo, apenas 16% dos afro-americanos foram classificados corretamente.
• Vamos tentar melhorar o modelo adicionando o endereço. O censo da Flórida tem os dados de endereço que serão adicionados ao modelo. Os dados estão em FLCensus.

• Variáveis:
  • county: id do distrito eleitoral no censo do condado
  • VTD: id do distrito eleitoral
  • total.pop: população do disitrito eleitoral
  • white: proporção de brancos não hispânicos no distrito eleitoral
  • black: proporção de afro-americanos não hispânicos no distrito eleitoral
  • api: proporção de asiáticos e das ilhas do Pacífico não hispânicos no distrito eleitoral
  • hispanic: proporção de hispânicos no distrito eleitoral
  • others: proporção de outros grupos raciais no distrito eleitoral

• Para calcular \(P(\mbox{race}|\mbox{surname e residence})\), nosso objetivo, vamos usar a seguinte versão da regra de Bayes: \[P(A|B\mbox{ e }C)=\frac{P(B|A\mbox{ e }C)P(A|C)}{P(B|C)}\]

• Aplicando a fórmula no nosso problema: \[P(\mbox{race}|\mbox{surname e residence})=\\\frac{P(\mbox{surname}|\mbox{race e residence})P(\mbox{race}|\mbox{residence})}{P(\mbox{surname}|\mbox{residence})}\]

• \(P(\mbox{race}|\mbox{residence})\) está disponível nos dados do censo da Flórida, as outras duas probabilidade condicionais, \(P(\mbox{surname}|\mbox{race e residence})\) e \(P(\mbox{surname}|\mbox{residence})\) não estão nos dados.
• Para encontrar essas probabilidade vamos fazer a hipótese de que o sobrenome e o endereço são independentes quando a raça é dada. Essa hipótese quer dizer que uma vez que sabemos a raça do eleitor, o endereço não nos dá informações adicionais a respeito do sobrenome.

• Com essa hipótese temos: \[P(\mbox{surname}|\mbox{race e residence})=\frac{P(\mbox{surname}|\mbox{race})P(\mbox{residence}|\mbox{race})}{P(\mbox{residence}|\mbox{race})}\\=P(\mbox{surname}|\mbox{race})\]

• Da lei da probabilidade total temos que: \[P(\mbox{surname}|\mbox{residence})=\\\sum_{\mbox{race}}P(\mbox{surname}|\mbox{race e residence})P(\mbox{race}|\mbox{residence})\\=\sum_{\mbox{race}}P(\mbox{surname}|\mbox{race})P(\mbox{race}|\mbox{residence})\]

• O censo tem dados para \(P(\mbox{race}|\mbox{surname})\), mas precisamos de conhecer \(P(\mbox{surname}|\mbox{name})\), felizmente já vimos a regra de Bayes: \[P(\mbox{surname}|\mbox{race})=\frac{P(\mbox{race}|\mbox{surname})P(\mbox{surname})}{P(\mbox{race})}\]

• Temos dados para \(P(\mbox{race}|\mbox{surname}\) e \(P(\mbox{surname})\), usaremos a lei da probabilidade total para obter \(P(\mbox{race})\): \[P(\mbox{race})=\sum_{\mbox{residence}}P(\mbox{race}|\mbox{residence})P(\mbox{residence})\]

• Vamos ao R! Começamos calculando \(P(\mbox{race})\), o primeiro passo é carregar os dados e ajustar a classificação das raças.

FLCensus <- read.csv("FLCensus.csv")

FLCensus <- FLCensus %>%
  rename("asian"="api", "other" = "others")

• Depois vamos usar a função weighted.mean() que calcula médias ponderadas.

race_prop <- FLCensus %>%
  select(total.pop, white, black, asian, hispanic, other) %>%
  pivot_longer(names_to = "race", values_to = "pct", -total.pop) %>%
  group_by(race) %>%
  summarize(race_prop = weighted.mean(pct, w=total.pop)) %>%
  arrange(desc(race_prop))

race_prop

## # A tibble: 5 × 2
##   race     race_prop
##   <chr>        <dbl>
## 1 white       0.579 
## 2 hispanic    0.225 
## 3 black       0.152 
## 4 asian       0.0242
## 5 other       0.0206

• Agora podemos calcular \(P(\mbox{surname}|\mbox{race})\):

FLVotersJoin <- left_join(FLVotersJoin, race_prop, by="race")
total.count <- sum(cnames$count)

cnames_reshape <- cnames %>%
  select(-race_pred, -highest_pct) %>%
  pivot_longer(names_to = "race", values_to = "pct", cols = starts_with("pct")) %>%
  mutate(race = str_replace(race, "pct", ""),
         race = recode(race, "api"="asian", "others"="other"),
         race_surname = pct/100) %>%
  select(-pct) %>%
  left_join(race_prop, by="race") %>%
  mutate(prob_surname_race = race_surname * (count/total.count)/race_prop) %>%
  rename("posib_race" = race)

• Conhecido \(P(\mbox{surname}|\mbox{race})\), podemos calcular \(P(\mbox{surname}|\mbox{residence})\)

merge_temp <- FLCensus %>%
  pivot_longer(names_to = "pop_race",
               values_to = "pop_pct",
               cols = c(white, black, hispanic, asian, other)) %>%
  inner_join(FLVoters, by=c("county", "VTD")) %>%
  inner_join(cnames_reshape, by=c("surname", "pop_race" = "posib_race")) %>%
  mutate(race_residence = prob_surname_race * pop_pct)

name_residence <- merge_temp %>%
  group_by(surname, county, VTD) %>%
  summarize(name_residence = sum(race_residence))

## `summarise()` has grouped output by 'surname', 'county'. You can override using
## the `.groups` argument.

• Agora podemos calcular \(P(\mbox{race}|\mbox{surname e residence})\)

FLVoters_full <- merge_temp %>%
  inner_join(name_residence, by=c("surname", "county", "VTD")) %>%
  mutate(pred_race = prob_surname_race * pop_pct/name_residence) %>%
  select(surname, pop_race, race, pred_race, county, VTD)

• Por fim, podemos comparar os resultados do modelo com a raça declarada. Assim como no caso anterior a raça prevista será a que teve maior probabilidade.

FL_updated <- FLVoters_full %>%
  ungroup() %>%
  group_by(surname, county, VTD) %>%
  filter(row_number(desc(pred_race))==1) %>%
  ungroup()

race_tp_new <- FL_updated %>%
  group_by(race) %>%
  summarize(tp = mean(race == pop_race)) %>%
  arrange(desc(tp))

## # A tibble: 6 × 2
##   race          tp
##   <chr>      <dbl>
## 1 white    0.942  
## 2 hispanic 0.857  
## 3 black    0.628  
## 4 asian    0.604  
## 5 other    0.00797
## 6 native   0

• A taxa de acerto, true positive, aumentou de 16% para 63% no grupo de afro-americanos e de 56% para 60% no grupo de asiáticos.
• A melhora nas previsões do modelo indicam que endereço é relevante para prever raça.

• Um exemplo pode ilustrar o efeito do endereço na previsão de raça. De acordo com os dados do censo 27% das pessoas com o sobrenome “White” são afro-americanos.

cnames %>%
  filter(surname == "WHITE") %>%
  pull()

## [1] 67.91

• Porém, ao considerarmos o endereço, essa probabilidade vai de praticamente 0 a 98%.

FLVoters_full %>%
  filter(surname == "WHITE", pop_race == "black") %>%
  select(pred_race) %>%
  summary()

##    pred_race       
##  Min.   :0.004588  
##  1st Qu.:0.072232  
##  Median :0.159496  
##  Mean   :0.250711  
##  3rd Qu.:0.293640  
##  Max.   :0.981864

• Para terminar o exemplo vamos checar os falsos positivos gerados pelo modelo.

race_fp_new <- FL_updated %>%
  group_by(pop_race) %>%
  summarize(fp = mean(race != pop_race)) %>%
  arrange(desc(fp))

race_fp_new

## # A tibble: 5 × 2
##   pop_race    fp
##   <chr>    <dbl>
## 1 other    0.778
## 2 asian    0.333
## 3 black    0.220
## 4 hispanic 0.212
## 5 white    0.122

• Os falsos positivos entre brancos foram de 20% para 12%. Grande parte dessa queda é devida à correção na classificação de afro-americanos que foram classificados como brancos no modelo que considera apenas o sobrenome.

• Uma variável aleatória atribui um valor numérico a cada evento em um experimento. Esses valores representam eventos mutuamente exclusivos e exaustivos, juntos formam todo o espaço amostral. Uma variável aleatória discreta toma um número finito, ou pelo menos contável, de valores distintos. Uma variável aleatória contínua toma um número de valores infinitos e não contáveis.

• Comecemos com o caso mais simples de variável aleatória discreta: jogar um moeda. Para esse experimento definimos uma variável aleatória binária, \(X\) que será igual a um se der cara e 0 caso contrário.
• Em geral, variáveis aleatórias que tomam dois valores disitintos são chamadas de variáveis aleatórias de Bernoulli. Isso vale para qualquer experimento com dois possíveis resultados distintos: \(\{\mbox{ganhar}, \mbox{perder}\}\), \(\{\mbox{subir}, \mbox{descer}\}\), \(\{\mbox{votar}, \mbox{não votar}\}\), etc.
• É comum dizermos que \(X=1\) em caso de sucesso e \(X=0\) em caso de falha e usarmos \(p\) para represetnar a probabilidade de sucesso.

• A distribuição de uma variável aleatória discreta pode ser caracterizada por uma função massa de probabilidade, ou simplesmente função de probabilidade, que será denotada por PMF (do inglês: probability mass function).
• A PMF, \(f(x)\), de uma variável aleatória \(X\) é definida como a probabilidade que a variável aleatória tome o valor \(x\), ou seja, \(f(x) = P(X=x)\).
• No caso de uma Bernoulli a PMF é dada por \(f(1)=P(X=1)=p\) e \(f(0)=P(x=0)=1-p\), para todos os outros valores \(f(x)=0\).

• Outra função importante é a função de distribuição acumulada, denotada por CDF (do inglês: cumulative distribution function). A CDF, \(F(x)\), representa a probabilidade acumulada dos valores iguais ou menores do que \(x\), ou seja, \(F(x) = P(x \leq X)\).
• A CDF representa a soma das PMF dos valores menores ou iguais a \(x\): \[F(x)=P(X \leq x) = \sum_{k \leq x}f(k)\]
• A CDF toma valores entre zero e um para qualquer variável aleatória, contínua ou discreta, e é uma função não decrescente porque não existe probabilidade negativa, de forma que quando \(x\) aumenta \(\sum_{k \leq x}f(k)\) não pode diminuir.

• A função de probabilidade de uma variável aleatória de **Bernoulli* com probabilidade de sucesso \(p\) é dada por: \[f(x)=\begin{cases}p \qquad \qquad \mbox{se x=1}\\1-p \ \qquad \mbox{se x=0}\\p \qquad \qquad \mbox{caso contrário}\end{cases}\]

• A função de densidade acumulada é dada por: \[F(x)=\begin{cases}0 \qquad \qquad \mbox{se x<0}\\1-p \ \qquad \mbox{se 0}\leq\mbox{x}<1\\1 \qquad \qquad \mbox{x}\geq1\end{cases}\]

• A variável aleatória uniforme é um exemplo simples de variável aleatória contínua. A variável aleatória uniforme atribui a mesma probabilidade para todos os número em um intervalo \([a, b]\), os números fora desse intervalo tem probabilidade zero.

• Como variáveis aleatórias contínuas podem tomar um número infinito e não contável de valores não podemos definir uma função de probabilidade, em vez disso usamos a função de densidade de probabilidade ou função de densidade, denotada por PDF (do inglês: probability density function), que quantifica a probabilidade que a variável aleatória \(X\) tome um valor \(x\).
• A PDF tem valores não negativos que podem ser maiores do que um, porém a área abaixo de uma PDF deve ser igual a um.

• A CDF de um variável aleatória contínua tem a mesma interpretação do que no caso de variáveis aleatórias discretas, \(F(x)=P(X\leq x)\), porém, no lugar da soma usamos a integral. \[F(x)=\int_{-\infty}^xf(t)dt\]
• Como a área abaixo da PDF é igual a um, temos que \(F(x)=1\) quando \(x \rightarrow \infty\).

• A função de densidade de uma variável aleatória uniforme no intervalo \([a,b]\) é dada por: \[f(x)=\begin{cases}\frac{1}{b-a} \qquad \ \mbox{se a}\leq x \leq \mbox{b}\\0 \qquad \ \ \ \ \ \mbox{caso contrário}\end{cases}\]

• A função de densidade acumulada é dada por: \[F(x)=\begin{cases}0 \qquad \qquad \ \mbox{se x} < \mbox{a}\\ \frac{x-a}{b-a} \qquad \quad \mbox{se a}\leq x < \mbox{b}\\1 \qquad \qquad \ \mbox{se x} \geq \mbox{b}\end{cases}\]

• A função dunif() calcula a PDF de uma distribuição uniforme (não é bem isso…).

dunif(0.5, min=0, max=1)

## [1] 1

dunif(1, min=-2, max=2)

## [1] 0.25

• A função punif() calcula a CDF de uma distribuição uniforme.

punif(1, min=-2, max=2)

## [1] 0.75

punif(0, min=-1, max=1)

## [1] 0.5

• Em uma distribuição uniforme entre zero e um a probabilidade de obter um valor menor do que \(x\) é igual a \(x\) para \(0 \leq x \leq 1\).

punif(0.6, min=0, max=1)

## [1] 0.6

punif(0.3, min=0, max=1)

## [1] 0.3

punif(0.8, min=0, max=1)

## [1] 0.8

• A função runif() retorna realizações de uma variável aleatória uniforme.

sims <- 1000
x <- runif(sims, min=0, max=1)
head(x, n=10)

##  [1] 0.01920359 0.40652059 0.56208492 0.95918199 0.77801306 0.07608611
##  [7] 0.33794735 0.74663095 0.53016555 0.78382456

• A função runif()e o resultado que em uma distribuição uniforme entre zero e um a probabilidade de obter um valor menor do que \(x\) é igual a \(x\) para \(0 \leq x \leq 1\) permitem simular uma variável aleatória Bernoulli.

p <- 0.5
y <- as.integer(x <= p)
head(y, n=10)

##  [1] 1 1 0 0 0 1 1 0 0 0

mean(y)

## [1] 0.507

• Simulação com \(p = 0,\!3\).

sims <- 1000
x <- runif(sims, min=0, max=1)
head(x, n=10)

##  [1] 0.2610260 0.6118536 0.9440752 0.9351696 0.4937785 0.3668355 0.8740632
##  [8] 0.9958626 0.5702721 0.5810730

p <- 0.3
y <- as.integer(x <= p)
head(y, n=10)

##  [1] 1 0 0 0 0 0 0 0 0 0

mean(y)

## [1] 0.31

• A distribuição Bernoulli é um caso particular da distribuição binomial. Em vez de jogar uma moeda uma vez, considere que a mesma moeda é jogada de forma independente várias vezes, a distribuição binomial representa o número de vezes que a moeda deu cara (sucessso).
• De forma mais geral, uma variável aleatória binomial, \(X\), representa o número de sucessos em \(n\) repetições independentes com probabilidade de sucesso igual a \(p\). A variável aleatória binomial pode tomar valores entre 0 e \(n\).

• Como a variável aleatória é discreta, a função de probabilidade representa a probabilidade de \(X\) tomar o valor \(x\), \(P(X=x)\).
• A CDF representa a probabilidade acumulada de uma variável aleatória binomial ter até \(x\) sucessos em \(n\) tentativas, \(P(X\leq x)\).

• A função de probabilidade de uma variável aleatória binomial com probabilidade de sucesso \(p\) e \(n\) repetições é dada por: \[f(x)=P(X=x)={n \choose x}p^x p^{n-x}\]

• A CDF é dada por: \[F(x)=P(X \leq x)= \sum_{k=0}^x{n \choose k}p^k p^{n-k}\]

• Por exemplo, em uma binomial com 3 repetições e probabilidade de sucesso igual a 0,5 a função de probabilidade para os valores possíveis de \(x\) será: \[f(0)={3 \choose 0}\times 0,\!5^0 \times (1-0,\!5)^{3-0}=\frac{3!}{(3-0)!0!}\times 0,\!5^3=0,\!125\] \[f(1)={3 \choose 1}\times 0,\!5^1 \times (1-0,\!5)^{3-1}=\frac{3!}{(3-1)!1!}\times 0,\!5^3=0,\!375\]

• Continuação… \[f(2)={3 \choose 2}\times 0,\!5^2 \times (1-0,\!5)^{3-2}=\frac{3!}{(3-2)!2!}\times 0,\!5^3=0,\!375\] \[f(3)={3 \choose 3}\times 0,\!5^3 \times (1-0,\!5)^{3-3}=\frac{3!}{(3-3)!3!}\times 0,\!5^3=0,\!125\]

• A CDF é dada por: \[P(X\leq 0)=F(0) = f(0) = 0,\!125\] \[P(X\leq 1) = F(1) = f(0)+ f(1) = 0,\!5\] \[P(X\leq 2) = F(2) = f(0)+ f(1)+f(2) = 0,\!875\] \[P(X\leq 3) = F(3) = f(0)+ f(1) +f(2)+f(3)= 1\]

• A função dbinom() retorna a PMF (probabilidade) de uma binomial

dbinom(2, size=3, prob=0.5)

## [1] 0.375

dbinom(0:3, size=3, prob=0.5)

## [1] 0.125 0.375 0.375 0.125

• A função pbinom() retorna a CDF de um binomial

pbinom(2, size=3, prob=0.5)

## [1] 0.875

pbinom(0:3, size=3, prob=0.5)

## [1] 0.125 0.500 0.875 1.000

• Outra distribuição importante, talvez a mais importante de todas, é a distribuição normal ou Gaussiana. Essa distribuição é particularmente importante porque para um grande número de variáveis aleatórias vale que a soma de muitas variáveis aleatórias com a mesma distribuição segue uma distribuição normal. Vamos tratar disso ainda nesta unidade.

• Uma variável aleatória normal pode tomar qualquer número real, \((-\infty, \infty)\). A distribuição normal é definida por dois parâmetros: média, \(\mu\), e desvio padrão, \(\sigma\).
• Se \(X\) é uma variável aleatória normal com média \(\mu\) e variância \(\sigma^2\)é comum escrevermos \(X \sim N(\mu, \sigma^2)\).

• A função de densidade, PDF, de uma variável aleatória normal com média \(\mu\) e variância \(\sigma^2\) é dada por: \[f(x)=\frac{1}{\sqrt{2\pi \sigma}}\exp\{-\frac{1}{2\sigma^2}(x-\mu)^2\}\]

• A CDF é dada por: \[F(x)=\int_{-\infty}^x f(t)dt= \int_{-\infty}^x\frac{1}{\sqrt{2\pi \sigma}}\exp\{-\frac{1}{2\sigma^2}(t-\mu)^2\}dt\]

• A PDF é simétrica em torno da média.
• Diferentes valores para a média deslocam a PDF sem mudar o formato.
• Valores mais altos para o desvio padrão deixam a PDF mais achatada e fazem com que a CDF cresça de forma mais gradual.

• Se \(X\) é uma distribuição normal com média 0 e desvio padrão igual a 1, \(x \sim N(0,1)\), diz-se que \(X\) é uma distribuição normal padrão.

• Suponha que \(X\) seja uma variável aleatória normal com média \(\mu\) e desvio padrão \(\sigma\), ou seja, \(X \sim N(\mu,\sigma^2)\). SEja \(c\) uma constante arbitrária. Vale que:
  • Uma variável aleatória definida por \(Z = X+c\) tem disitribuição normal, \(Z \sim N(\mu+c, \sigma^2)\).
  • Uma variável aleatória definida por \(Z = cX\) tem distribuição normal, \(Z \sim N(c\mu, (c\sigma)^2\).

• Essas propriedade implicam que o z-score de uma variável aleatória normal segue uma distribução normal padrão: \[\mbox{z-score}=\frac{X-\mu}{\sigma}\sim N(0,1)\]

• Se os dados tem distribuição normal padrão, então cerca de 2/3 dos dados estão entre a média menos um desvio padrão e a média mais um desvio padrão. Da mesma forma, cerca de 95% dos dados estão entre dois desvios padrão da média.
• Como desvio padrão é 1 e a média é 0, cerca de 2/3 dos dos dados estão entre -1 e 1 e cerca de 95% dos dados estão entre -2 e 2.

• Para confirmar essa propriedade vamos encontrar a probabilidade de uma variável aleatória com média \(\mu\) e desvio padrão \(\sigma^2\) estar a \(k\) desvios padrão da média, \(k>0\). Para simplificar vamos considerar o z-score que tem uma distribuição normal padrão. \[P(\mu - k\sigma \leq X \leq \mu + k\sigma)=P(-k\sigma \leq X-\mu \leq k\sigma)\\=P(-k \leq \frac{X-\mu}{\sigma} \leq k)=P(-k \leq Z \leq k)\]
• Onde \(Z\sim N(0,1)\).

• Agora, repare que:

\[P(-k \leq Z \leq k)=P(Z \leq k) - P(Z \leq -k)= F(k)-F(-k)\]

• A próxima figura ilustra o resultado. A área verde, entre \(-k\) e \(k\), é a que desejamos, ou seja, \(P(-k \leq Z \leq k)\).
  • A área cinza, a esquerda de \(-k\), é \(P(Z \leq -k)\).
  • A soma da área verde com a área cinza é \(P(Z \leq k)\).
  • Logo, a área deseja é a soma da área verde com a área cinza menos a área cinza, ou seja, \(P(Z \leq k) - P(Z \leq -k)\).

• Podemos usar o R para confirmar que em uma normal padrão cerca de 2/3 da área abaixo da PDF está entre -1 e 1, ou a um desvio padrão da média, e cerca de 95% da área abaixo da PDF está entre -2 e 2, ou a dosi desvios padrão da média.

pnorm(1)-pnorm(-1)

## [1] 0.6826895

pnorm(2)-pnorm(-2)

## [1] 0.9544997

• Usamos normal padrão para simplificar o argumento, mas o resultado vale para qualquer normal.

mu <- 5
sigma <- 2

pnorm(mu+sigma, mean=mu, sd = sigma) - pnorm(mu-sigma, mean=mu, sd=sigma)

## [1] 0.6826895

pnorm(mu+2*sigma, mean=mu, sd = sigma) - pnorm(mu-2*sigma, mean=mu, sd=sigma)

## [1] 0.9544997

• Como exemplo de uso da distribuição normal considere novamente o fênomeno da regressão para média que discutimos na unidade 4.
• Lá apresentamos evidências de que em estados que Obama recebeu uma grande fração de votos em 2008 era provável que ele recebesse uma fração menor em 2012.
• Para ilustrar o ponto fizemos uma regressão entre os votos de Obama e 2008 e 2012.

pres08 <- read.csv("pres08.csv")
pres12 <- read.csv("pres12.csv")

pres <- full_join(pres08, pres12, by="state")

pres <- pres %>%
  mutate(Obama2008.z = as.numeric(scale(Obama.x)),
         Obama2012.z = as.numeric(scale(Obama.y)))

fit1 <- lm(Obama2012.z ~ -1+Obama2008.z, data=pres)
fit1

## 
## Call:
## lm(formula = Obama2012.z ~ -1 + Obama2008.z, data = pres)
## 
## Coefficients:
## Obama2008.z  
##      0.9834

• Vamos examinar os resíduos da regressão em fit1 e comparar com uma distribuição normal. Primeiros faremos o histrograma dos resíduos e a densidade da normal, depois usaremos o QQ-plot.
• Para que os resíduos sejam comparáveis com uma distribuição normal padrão usaremos a função scale() nos resíduos. De forma alternativa podemos dividir os resíduos pelo desvio padrão dos resíduos (lembre que os resíduos tem média zero).
• Os gráficos serão feitos com R básico.

e <- resid(fit1)
e.zscore <- scale(e)

hist(z.score, freq=FALSE, ylim=c(0,4),
     xlab = "Standardized residuals",
     main = "Distribution of standardized residuals")
x <- seq(from=-3, to=3, by=0.01)
lines(x, dnorm(x))

qqnorm(e.zscore, xlim=c(-3,3), ylim=c(-3,3))
abline(0,1)

• Agora considere o seguinte modelo: \[\mbox{Votos padronizados para Obama em 2012}=\\0,\!98\times\mbox{Votos padronizados para Obama em 2008} + \varepsilon\]
• Onde \(\varepsilon \sim N(0, 0,\!18)\).

• Lembre que a média de \(\varepsilon\) é zero por construção, podemos usar o R para encontrar o desvio padrão e confirmar que a média é zero.

mean(e)

## [1] -1.093832e-16

e.sd <- sd(e)
e.sd

## [1] 0.1812239

• O modelo proposto descreve um potencial processo gerador de dados para a fração de votos recebidos por Obama em 2012.
• Considere a Califórnia onde Obama teve 61% dos votos em 2008, o valor padronizado é de 0,87.

(61-mean(pres$Obama.x))/sd(pres$Obama.x)

## [1] 0.8720631

• De acordo com o modelo, qual a probabilidade de Obama ter mais de 61% dos votos da Califórnia em 2012?
• A função dnorm() responde nossa pergunta. Como a função retorna a probabilidade da variável aleatória ser menor do que \(x\) e queremos a probabilidade de ser maior vamos calcular 1 menos a probabilidade de ser menor que \(x\), ou, mais fácil, usar o argumento lower.tail como falso na função.

CA.2008 <- filter(pres, state == "CA") %>%
  select(Obama2008.z) %>%
  pull()

CA.mean2012 <- coef(fit1) * CA.2008
CA.mean2012

## Obama2008.z 
##   0.8576233

1 - pnorm(CA.2008, mean=CA.mean2012, sd=e.sd)

## [1] 0.4682463

pnorm(CA.2008, mean=CA.mean2012, sd=e.sd, lower.tail = FALSE)

## [1] 0.4682463

• A probabilidade de em 2012 Obama conseguir uma proporção dos votos maior do que em 2008 é de 46.8%. Façamos o mesmo exercíico para o Texas onde Obama recebeu 44% dos votos, ou -0,67 se considerarmos os valores padronizados.

TX.2008 <- filter(pres, state == "TX") %>%
  select(Obama2008.z) %>%
  pull()

TX.mean2012 <- coef(fit1) * TX.2008
TX.mean2012

## Obama2008.z 
##  -0.6567543

pnorm(TX.2008, mean=TX.mean2012, sd=e.sd, lower.tail = FALSE)

## [1] 0.5243271

• No estado onde Obama foi bem em 2008, a probabilidade de melhorar em 2012, 46,8%, é menor do que a probabilidade de piorar. Já no estado onde Obama foi mal em 2008, a probabilidade de melhorar em 2012, 52,4%, é maior do que a de piorar.
• Essa é a ideia de regressão à média.

• Funções como a função de probabilidade, função de densidade e função de distribuição acumulada descrevem completamente a disitribuição de uma variável aleatória, porém algumas vezes precisamos de informações mais compactas sobre a distribuição de uma variável aleatória.
• Uma medida útil nessas ocasiões é o valor esperado. O valor esperado pode ser visto como a média de uma variável aleatória, mas não deve ser confundido com a média amostral. Essa última diz respeito ao valor médio de uma amostra, enquanto o valor esperado represetna o valor médio sob um determinada distribuição.
• A média amostral muda de uma amostra para outra, o valor esperado é um conceito teórico e é fixo para uma dada distribuição.

• Considere, por exemplo, uma distribuição Bernoulli com probabilidade de sucesso igual a \(p\). Qual o valor esperado desta distribuição?
• Como ela só toma dois valores, 0 ou 1, o valor esperado é calculado como a média ponderada desses valores com pesos \(p\) e \(1-p\). \[E(X)=0\times P(X=0)+1\times P(X=1)=0\times f(0)+1\times f(1)\\=0\times (1-p) + 1\times p=p\]

• No caso de uma binomial, que pode tomar qualquer inteiro não negativo menor ou igual a \(n\), ou seja, pode tomar valores em \(\{0,1,\dots,n\}\), o valor esperado é dado por: \[E(X)=0\times f(0) + 1\times f(1)+\cdots+n\times f(n)=\sum_{x=0}^n x\times f(x)\]

• A ideia de soma ponderada vale para variáveis contínuas, porém, como essas variáves tomam um número infinito e não contável de valores, será preciso substituir a soma pela integral.
• Para uma uniforme no intervalo \([a b]\) o valor esperado é dado por: \[E(X)=\int_a^b x\times f(x)dx=\int_a^b \frac{x}{b-a}dx=\frac{x^2}{2(b-a)}\Bigg|_a^b=\frac{a+b}{2}\]

• O valor esperado ou esperança de uma variável aleatória é denotado por \(E(X)\) e definido por: \[E(X)=\begin{cases}\sum_x x\times f(x) \qquad \qquad \mbox{se }X\mbox{ for discreta}\\\int x\times f(x)dx\qquad \quad \ \ \ \mbox{se }X\mbox{ for contínua}\end{cases}\]
• Onde \(f(x)\) é afunção de probabilidade (função de densidade) de uma variável aleatória discreta (contínua).

• Sejam \(X\) e \(Y\) variáveis aleatória e \(a\) e \(b\) constantes arbitrárias. O valor esperado é um operador linear que satisfaz as seguintes propriedades:
  • \(E(a)=a\)
  • \(E(a X)=a E(X)\)
  • \(E(aX + b) = a E(X) + b\)
  • \(E(aX + bY) = a E(X) + bE(Y)\)
  • Se \(X\) e \(Y\) são independentes, então \(E(XY)=E(X)\times E(Y)\), mas, em geral, \(E(XY) \neq E(X) \times E(Y)\).

• Podemos usar as propriedades acima para encontrar o valor esperado de uma binomial de tamanho \(n\) e probabilidade de sucesso igual a \(p\).
• Para isso perceba que um binomial deste tipo é dada pela soma de \(n\) Bernoullis independentes com probabilidade de sucesso igual a \(p\), logo, chamando a binomial de \(X\) e as Bernoulli de \(Y_i\): \[E(X)= E\left(\sum_{i=1}^n Y_i \right)=\sum_{i=1}^nE(Y_i)=np\]

• Outra estatística importante é o desvio padrão (e a variância, que é o quadrado do desvio padrão) de uma variável aleatória. Assim como no valor esperado, é importante diferenciar o desvio padrão de uma amostra do conceito teórico de desvio padrão de uma variável aleatória.

• A variância de uma variável aleatória \(X\) é definida como: \[V(X)=E\{[X-E(X)]^2\}\]
• O desvio padrão é a raiz quadrada da variância.

• Usando as propriedades do valor esperado podemos mostrar que a variância é a diferença entre a esperança de \(X^2\) e o quadrado da esperança de \(X\): \[V(X)=E\{[X-E(X)]^2\}=E\{X^2 - 2XE(X)+[E(X)]^2\}\\=E(X^2)-2E(X)E(X)+[E(X)]^2\\=E(X^2)-[E(X)]^2\]
• A esperança de \(X^2\) é chamada de segundo momento de \(X\) e a esperança de \(X\) é chamada de primeiro momento.

• A expressão acima pode ser usada para calcular a variância de uma Bernoulli com probabilidade de sucesso dada por \(p\). Para isso note que se \(X\) tem essa distribuição, então \(X^2=X\), pois \(0^2=0\) e \(1^2=1\). \[V(X)=E(X^2)-[E(X)]^2=p-p^2=p(1-p)\]

• Também podemos calcular a variância de uma uniforme no intervalo \([a b]\).

\[V(X)=E(X^2)-[E(X)]^2=\int_a^b \frac{x^2}{b-a} - \left(\frac{a+b}{2} \right)^2\\=\frac{x^3}{3(b-a)}\Bigg|_a^b-\left(\frac{a+b}{2} \right)^2=\frac{1}{12}(b-a)^2\]

• As fórmulas acima podem ser usadas para comparar a variância de uma variável aleatória com a variância de uma amostra obtida com realizações desta variável aleatória.

• Comecemos com a Bernoulli.

p <- 0.5

p * (1-p)

## [1] 0.25

y <- runif(500, 0,1) >= p
var(y)

## [1] 0.2502445

• Agora com uma uniforme no intervalo \([1 5]\)

a <- 1
b <- 5

(1/12)*(b-a)^2

## [1] 1.333333

var(runif(500, min=1, max=5))

## [1] 1.321897

• Sejam \(X\) e \(Y\) variáveis aleatórias e \(a\) e \(b\) constantes arbitrárias. A variância é um operador com as seguintes propriedades:
  • \(V(a)=0\)
  • \(V(aX)=a^2V(X)\)
  • \(V(X+b)=V(X)\)
  • \(V(aX+b)=a^2 V(x)\)
  • Se \(X\) e \(Y\) são independentes \(V(X+Y)=V(X)+V(Y)\).

• Com as propriedades acima podemos calcular a variância de uma binomial com tamanho \(n\) e probabilidade de sucesso igual a \(p\). \[V(X)=V\left(\sum_{i=1}^n Y_i\right)=\sum_{i=1}^nV(Y_i)=np(1-p)\]

• Outro exemplo é a soma de duas normais independentes (lembre que a soma de normais é uma normal), especificamente: seja \(X\sim N(\mu_X, \sigma_X^2)\), \(Y\sim N(\mu_Y, \sigma_Y^2)\) e \(Z=aX + bY + c\): \[E(Z)=E(aX + bY + c)=aE(X)+bE(Y)+c=a\mu_X + b\mu_Y + c\\V(Z)=V(aX + bY)=a^2V(X)+b^2V(Y)=a^2 \sigma_X^2 + b^2 \sigma_Y^2\]
• Logo, \(Z \sim N(a\mu_X + b\mu_Y + c,a^2 \sigma_X^2 + b^2 \sigma_Y^2)\)

• Na unidade 4 discutimos um modelo de previsão para eleições. O problema é que quando trabalhamos com amostras estamos limitados às estatísticas observadas na amostra(e.g. média amostral) e essas mudam de amostra para amostra, isso se chama variabilidade amostral.

• Em aplicações não conhecemos os parâmetros da distribuição da variável aleatória que gerou os dados, apenas podemos calcular as versões amostrais, porém, essas mudam de amostra para amostra.
• Como saber se estamos certos, ou pelo menos quão perto de estarmos certos? Vamos voltar ao modelo de previsão de resultados eleitorais para discutir essa questão.

• Suponha que a fração de votos para candidatos dos partidos que não o Republicano e Democrata é tão pequena que pode ser ignorada, ou seja vamos tratar a eleição com apenas dois resultados possíveis.
• Suponha também que a “verdadeira” proporção de votos para Obama e McCain será revelada na apuração de votos. Essa proporção muda de estado para estado, de forma que no estado \(j\) a proproção de Obama é \(p_j\) e a de McCain é \(1-p_j\).

• Os dados estão em pres08, vamos criar \(p\) dividindo os votos de Obama pela soma de votos de Obama e McCain.

pres08 <- read.csv("pres08.csv")

pres08 <- pres08 %>%
  mutate(p = Obama/(Obama + McCain))

• Vamos criar amostras hipóteticas de 1000 eleitores (pense em pesquisas com 1000 eleitores). Para cada estado temos uma distribuição binomial com probbailidade de sucesso \(p_j\) e tamanho 1000.
• Usando Monte Carlo vamos estimar a proporção de votos para Obama em cada estado, \(p_j\). Se \(p_j\) for maior do que 0,5 os votos do estado no colégio eleitoral vão para Obama, se for menor vão para McCain.
• Por conta da variabilidade amostral repetiremos o exercício muitas vezes.

• Para fazer as simulações da binomial vamos usar a função rbinom() com \(n=1000\) e \(p_j\) dado por um vetor com as probabilidades de cada estado.

• A função vai retornar o número de votos em 1000 associado a cada elemento do vetor.

rbinom(2, 1000, prob = c(0.1, 0.9))

## [1]  80 893

• Vamos criar uma função que terá como argumento o número de repetições do exercício, o data.frame com os dados e número de eleitores em cada amostra e retorna os votos do Obama no colégio eleitoral.

sim_election <- function(.id, df, n_draws=1000){
  mutate(df, draws = rbinom(n(), n_draws, p)) %>%
    filter(draws > n_draws/2) %>%
    summarize(EV=sum(EV), .id=.id)
}

sim_election(1, pres08, n_draws = 1000)

##    EV .id
## 1 364   1

sim_election(2, pres08, n_draws = 1000)

##    EV .id
## 1 353   2

• Para repetir os exercícios usaremos a função map_df() do pacote purr que é parte do tidyverse.

sims <- 10000
sim_results <- map_df(seq_len(sims), ~ sim_election(.x, pres08, n_draws = 1000))

• O data.frame sim_results tem o número de votos para Obama no colégio eleitoral a cada repetição do exercício.

head(sim_results)

##    EV .id
## 1 326   1
## 2 375   2
## 3 367   3
## 4 364   4
## 5 364   5
## 6 338   6

• Vamos usar um gráfico para ilustrar o resultado das simulações. O gráfico vai mostrar o histograma de todas as simulações e marcar com uma linha azul o número de votos que Obama de fato teve no colégio eleitoral.

sim_results %>%
  ggplot(aes(EV, y= after_stat(density))) +
  geom_histogram(binwidth = 12) +
  geom_vline(xintercept = 364, color = "blue", linewidth=0.5) +
  labs(x="Obama's Electoral College votes", y="Density") +
  theme_classic()

• O número de votos no colégio eleitoral muda em cada exercício, isso ocorre por conta da variabilidade amostral, mas em todos os exercícios Obama ficou com mais de 270 votos, quantidade necessária para ser eleito.
• Vamos dar uma olhada em algumas estatísticas dos resultados da simulação.

sim_results %>%
  select(EV) %>%
  summarise_all(list(mean = mean,
                     median = median,
                     var = var,
                     sd = sd,
                     max = max,
                     min=min))

##       mean median      var       sd max min
## 1 352.3853    353 270.0883 16.43436 393 285

• Podemos calcular analiticamente o valor esperado do votos para Obama no Colégio Eleitoral gerado por esse modelo.

• Seja \(S_j\) o número de entrevistados, em um total de mil, que responderam votar em Obama no estado \(j\) e \(v_j\) o número de votos do estado no Colégio Eleitoral. \[E(\mbox{Votos para Obama})=\sum_{j=1}^{51} v_j \times P(\mbox{Obama vence no estado } j)\\=\sum_{j=1}^{51} v_j\times P(S_j > 500)\]

• Para calcular as probabilidade usaremos a função pbinom(), como queremos maior ou igual uaremos o argumento lower.tail como falso.

n_draws <- 1000

pres08 %>%
  mutate(pb = pbinom(n_draws/2, size=n_draws, prob=p, lower.tail=FALSE)) %>%
  summarise(mean = sum(pb*EV))

##       mean
## 1 352.1388

• Também podemos calcular analiticamente a variância e o desvio padrão. \[V(\mbox{Votos para Obama})=\sum_{j=1}^{51} v_j^2 P(S_j > 500)(1-P(S_j > 500))\]

pres08 %>%
  mutate(pb = pbinom(n_draws/2, size=n_draws, prob=p, lower.tail=FALSE)) %>%
  summarize(V = sum(pb*(1-pb)*EV^2),
            sd = sqrt(V))

##          V       sd
## 1 268.8008 16.39515

• Suponha uma amostra aleatória de \(n\) independente, identicamente distribuídas observações, \(X_1, X_2, \dots, X_n\) de uma disitribuição de probabilidade com valor esperado \(E(X)\). A lei dos grandes números é dada por: \[\bar{X}_n = \frac{1}{n}\sum_{i=1}^n X_i \rightarrow E(X)\]
• A lei dos grandes números diz que a medida que o tamanho da amostra aumenta, a média amostral converge para o valor esperado.

• Para ilustrar a lei dos grande números vamos fazer uma simulação de Monte Carlo. Construiremos amostras aleatórias de uma distribuição binomial com probabilidade de sucesso igual a 0,2 e tamanho 10. Sabemos que o valor esperado é 2 (\(E(X)=np=10\times 0,\!2=2\)), vamos calcular a média amostral para amostras de diversos tamanhos.
• Quanto maior o tamanho da amostra, mais perto de 2 deve ser a média amostral.

• Para a simulação começamos com uma amostra de tamanho dois, depois vamos tirando mais números e com isso aumentando o tamanho da amostra. Repetimos o processo até termos uma amostra de 1000 observações.

• Para calcular as médias usaremos a função cumsum() que calcula a soma cumulativa de um vetor.

cumsum(c(1,2,3,4))

## [1]  1  3  6 10

cumsum(1:10)

##  [1]  1  3  6 10 15 21 28 36 45 55

sims <- 5000
p <- 0.2
size <- 10

lln_bin <- tibble(n = seq_len(sims),
                  x = rbinom(sims, prob=p, size=size),
                  mean = cumsum(x)/n,
                  distrib = str_c("Binomial (", size, ", ", p, ")"))

head(lln_bin, n=10)

## # A tibble: 10 × 4
##        n     x  mean distrib           
##    <int> <int> <dbl> <chr>             
##  1     1     1  1    Binomial (10, 0.2)
##  2     2     3  2    Binomial (10, 0.2)
##  3     3     2  2    Binomial (10, 0.2)
##  4     4     4  2.5  Binomial (10, 0.2)
##  5     5     4  2.8  Binomial (10, 0.2)
##  6     6     3  2.83 Binomial (10, 0.2)
##  7     7     1  2.57 Binomial (10, 0.2)
##  8     8     1  2.38 Binomial (10, 0.2)
##  9     9     2  2.33 Binomial (10, 0.2)
## 10    10     2  2.3  Binomial (10, 0.2)

tail(lln_bin, n=10)

## # A tibble: 10 × 4
##        n     x  mean distrib           
##    <int> <int> <dbl> <chr>             
##  1  4991     2  1.99 Binomial (10, 0.2)
##  2  4992     4  1.99 Binomial (10, 0.2)
##  3  4993     2  1.99 Binomial (10, 0.2)
##  4  4994     4  1.99 Binomial (10, 0.2)
##  5  4995     2  1.99 Binomial (10, 0.2)
##  6  4996     3  1.99 Binomial (10, 0.2)
##  7  4997     2  1.99 Binomial (10, 0.2)
##  8  4998     3  1.99 Binomial (10, 0.2)
##  9  4999     3  1.99 Binomial (10, 0.2)
## 10  5000     1  1.99 Binomial (10, 0.2)

• Vamos fazer um gráfico para ilustrar o processo de convergência.

lln_bin %>%
  ggplot() +
  geom_line(aes(x=n, y=mean)) +
  geom_hline(yintercept = 2, lty= "dashed") +
  labs(title = "Binomial (10, 0,2)", x="Sample size", y="Sample mean") +
  theme_classic()

• Podemos fazer exercício semelhante para uma distribuição uniforme [0, 1], nesse caso o valor esperado é 0,5.

lln_unif <- tibble(n = seq_len(sims),
                  x = runif(sims),
                  mean = cumsum(x)/n,
                  distrib = str_c("Uniforme (0, 1)"))

head(lln_unif, n=10)

## # A tibble: 10 × 4
##        n     x  mean distrib        
##    <int> <dbl> <dbl> <chr>          
##  1     1 0.679 0.679 Uniforme (0, 1)
##  2     2 0.143 0.411 Uniforme (0, 1)
##  3     3 0.166 0.329 Uniforme (0, 1)
##  4     4 0.884 0.468 Uniforme (0, 1)
##  5     5 0.557 0.486 Uniforme (0, 1)
##  6     6 0.220 0.441 Uniforme (0, 1)
##  7     7 0.358 0.430 Uniforme (0, 1)
##  8     8 0.230 0.405 Uniforme (0, 1)
##  9     9 0.838 0.453 Uniforme (0, 1)
## 10    10 0.538 0.461 Uniforme (0, 1)

tail(lln_unif, n=10)

## # A tibble: 10 × 4
##        n      x  mean distrib        
##    <int>  <dbl> <dbl> <chr>          
##  1  4991 0.544  0.499 Uniforme (0, 1)
##  2  4992 0.847  0.499 Uniforme (0, 1)
##  3  4993 0.338  0.499 Uniforme (0, 1)
##  4  4994 0.272  0.499 Uniforme (0, 1)
##  5  4995 0.0438 0.499 Uniforme (0, 1)
##  6  4996 0.853  0.499 Uniforme (0, 1)
##  7  4997 0.890  0.499 Uniforme (0, 1)
##  8  4998 0.335  0.499 Uniforme (0, 1)
##  9  4999 0.0394 0.499 Uniforme (0, 1)
## 10  5000 0.268  0.499 Uniforme (0, 1)

• Mais uma vez um gráfico pode ilustrar o processo de convergência.

lln_unif %>%
  ggplot() +
  geom_line(aes(x=n, y=mean)) +
  geom_hline(yintercept = 0.5, lty= "dashed") +
  labs(title = "Uniforme (0, 1)", x="Sample size", y="Sample mean") +
  theme_classic()

• Suponha que temos uma amostra aleatória de \(n\) i.i.d observações \(X_1, X_2, \dots, X_n\) de uma distribuição de probabilidade com valor esperado \(E(X)\) e variância \(V(X)\). Vamos denotar a média amostral por \(\bar{X}_n = \frac{\sum_{i=1}^n X_i}{n}\). O teorema central do limite diz que: \[\frac{\bar{X}_n - E(X)}{\sqrt{V(X)/n}} \leadsto N(0,1)\]
• O símbolo \(\leadsto\) indica convergência em distribuição a medida que \(n\) cresce.

• O teorema central do limite diz que o z-score da média amostral converge em distribuição para uma normal padrão a medida que o tamanho da amostr aumenta.
• Esse teorema é fundamental porque, mesmo que não conheçamos a distribuição verdadeira que gerou nossos dados, sabemos que para um amplo conjunto de distribuição a média amostral converge para uma normal.

• Assim como fizemos com a lei dos grandes números, vamos ilustrar o teorema central do limite com simulações de Monte Carlo.

sims <- 1000
n.samp <- 1000
z.binom <- z.unif <- rep(NA, sims)

for (i in 1:sims){
  x <- rbinom(n.samp, p=0.2, size=10)
  z.binom[i] <- (mean(x)-2)/sqrt(1.6/n.samp)
  x <- runif(n.samp, min=0, max=1)
  z.unif[i] <- (mean(x)-0.5)/sqrt(1/(12*n.samp))
}

results <- tibble(z.binom = z.binom,
                  z.unif = z.unif,
                  n.samp = seq(1:n.samp))

• Para melhor observar os resultados faremos um gráfico com o histograma da amostra padronizada e com a densidade de uma N(0, 1).
• Cada gráfico será feito com dois tamanhos de amostra.

• Gráfico para Binomial

results %>%
  ggplot() +
  geom_histogram(aes(x=z.binom, y=after_stat(density)), bins=20) +
  stat_function(fun = dnorm, color = "blue") +
  labs(title = "Binomial (0,2, 10)",
       x="z-score", y="Density") +
  theme_classic()

• Binomial, amostra com 100

• Binomial, amostra com 1000

• Gráfico para Uniforme

results %>%
  ggplot() +
  geom_histogram(aes(x=z.unif, y=after_stat(density)), bins=20) +
  stat_function(fun = dnorm, color = "blue") +
  labs(title = "Uniforme (0, 1)",
       x="z-score", y="Density") +
  theme_classic()

• Uniforme, amostra com 100

• Uniforme, amostra com 1000