A Moeda Viciada

Observação: Favor encaminharem para o email henriquealvarenga@gmail.com quaisquer erros ou sugestões no script ou na análise a seguir. Muito obrigado pelas correções e sugestões.

Objetivos de Aprendizagem

Objetivos Gerais

  • Compreender conceitos de probabilidade
  • Compreender o uso da curva normal teórica no cálculo do valor de p

Objetivos Específicos

  • Compreender o uso da distribuição nula no cálculo do valor de p
  • Compreender o que é uma distribuição amostral
  • Compreender o conceito de amostragem probabilística
  • Saber fazer simulações simples no r com a função sample()
  • Saber calcular o valor de p usando a função pnorm()

O Problema da Moeda Viciada

Você apostou com um colega quem iria tirar mais coroas num jogo de moedas e acabou perdendo. Depois de jogar a moeda 1000 vezes o resutado foi e 455 coroas e 545 caras. Ou seja, o nº de caras (aposta do colega) foi maior que o nº de coroas (sua aposta). Entretanto, você ficou intrigado com esse resultado e resolver verificar qual a probabilidade desse resultado ter sido encontrado. Chegando em casa pesquisou na internet como fazer esse cálculo usando o RStudio. Sua pesquisa foi toda documentada abaixo, pois você pretende confrontar seu colega com esses dados.

Formulação das hipóteses

  • H0 - Hipótese Nula: a moeda é honesta, ou seja, nº caras = 500
  • H1 - Hipótese Alternativa: a moeda é viciada, ou seja, nº caras > 500

Parte 1 - Aprendendo a fazer uma simulação no R

Vamos simular um jogo de moedas no R. A primeira coisa a fazer é criar um objeto com as propriedades que nos interessam de uma moeda, ou seja, possuir cara e coroa. veja como isso é simples no R. O comando novo aqui é o c(), que significa “combine”, combine. Esse comando combina os elementos dentro do parenteses.

moeda <- c("cara", "coroa")

O que temos de fazer agora é jogar a moeda. Para isso usamos a função sample(), que tem como argumentos o objeto moeda e o número de itens que iremos amostrar (escolher) ao jogar a moeda, isto é, o número de amostras que iremos retirar desse conjunto moeda.

Nesse exemplo, iremos retirar apenas 1 amostra do conjunto de 2 possibilidades. Em outras palavras: temos 2 possibilidades, cara ou coroa, ao simular uma jogada da moeda teremos UM unico resultado. A simulação, portanto, é de retirar ao acaso UMA das duas possibilidades.

sample(moeda, size=1, replace=TRUE)
## [1] "coroa"

Se desejarmos repetir essa jogada milhares de vezes, podemos usar o mesmo comando sample, mas definindo o parâmetro REPLACE com sendo TRUE. Esse parâmetro determina que após retirar um elemento da amostra, esse elemento volte para o conjunto antes que seja retirada a próxima amostra. Vamos colocar o resultado dessas jogadas num outro objeto, que chamaremos de mil.jogadas Verifique digitando o nome do objeto: mil.jogadas.

mil.jogadas <- sample(moeda, size=1000, replace = TRUE)
head(mil.jogadas)
## [1] "coroa" "coroa" "cara"  "cara"  "cara"  "cara"

Vamos verificar esses dados com um gráfico de barras:

barplot(table(mil.jogadas), col="lightblue")

Parte 2 - Fazendo a simulação para definir se a moeda é viciada

Vamos fazer o seguinte, vamos definir uma nova moeda, no qual cara=1, e coroa =0

moeda <- c(1, 0) # em nossa moeda,  1=cara, 0=coroa

Vamos gora jogar a moedas milhares de vezes e contar a quantidade de vezes em que o resultado foi “cara”.

sum(sample(moeda, size=1000, replace = TRUE)) # joga a moeda 1000 vezes, soma o resultado
## [1] 505

A Distribuição Amostral

A grande dificuldade da análise estatística é compreender a necessidade do uso de uma distribuição amostral. Para solucionar essa dificuldade é preciso compreender que a aposta que foi feita, na qual foram feitas 1000 jogadas e calculado a SOMA das vezes em que o resultado foi cara, é apenas UMA AMOSTRA dentro de um universos de infinitas apostas.

Precisamos imaginar um cenário fictício no qual infinitas apostas foram feitas, sendo que em cada aposta foram jogadas 1000 vezes as moedas e calculadas quantas caras foram obtidas em cada aposta.

A questão a ser respondida é justamente a seguinte: nesse universos de infinitas apostas com 1000 jogadas, qual a possibilidade de que tenha sido obtido esse resultado de 545 caras.

Esse universo de 1000 apostas é nossa distribuição amostral. Se a moeda for honesta, a média de caras deveria, teoricamente, ser 500 e não 545.

E então, vamos repetir por dez mil vezes essas mil jogadas e a contagem das caras em cada uma. Teremos ao final um conjunto de dez mil números que representam, cada um, o número de caras em cada uma das mil jogadas.

# recriando a moeda (essa linha seria denecessária, pois já criamos a moeda, mas preferi deixar aqui assim mesmo)
moeda <- c(1, 0) # em nossa moeda,  1=cara, 0=coroa
# a função sample com size =1000 joga a moeda 1000 vezes
# a função sum soma a quantidade de caras ( de nº "1")
# a função replicate(10000, ...) repete essa operação dez mil vezes
n.caras <- replicate(10000,sum(sample(moeda, size=1000, replace = TRUE)))

Vamos agora construir um histograma da distribuição amostral da soma do número de caras (n.caras):

hist(n.caras, col="lightblue", xlab ="nº de vezes que o resultado foi cara")

O valor esperado da média desses dados deve ser 500. Veja bem, se jogamos a moeda 1.000 vezes, teoricamente deveriam haver 500 caras e 500 coroas. Certamente que nem sempre será assim.

Entretanto, mas depois de repetir essa operação dez mil vezes é esperado que a média esteja bem próximo do valor de 500. O histograma nos mostra que a média está realmente ao redor de 500. Mas vamos conferir essa inspeção visual calculando a média e desses dados com a função mean() e vamos colocar esses resultados na variável m.caras:

m.caras   <- mean(n.caras)
m.caras
## [1] 500.1612

Calculando o nº de vezes em que houve mais de 545 caras em cada uma das 1000 jogadas.

up <- sum(n.caras >= 545)
up
## [1] 26

Calculando a proporção de vezes em que houve mais de 545 caras

up/10000
## [1] 0.0026

Essa proporção é a probabilidade de que em mil jogadas da moeda você obtenha o resultado caras 545 vezes ou mais. Como essa proporção foi muito pequena, ou seja, foi menor que o valor habitualmente usado de 0.05, podemos questionar nosso colega alegando que a meda estava provavelnmente viciada, pois é muito pouco provável obter esse resultado numa moeda honesta.

Calculando o Valor de P com a curva normal teórica

O cálculo da probabilidade de haver um resultado igual ou maior a 545 pode ser feito usando a curva normal teórica. Para isso precisamos saber apenas os parâmetros da distribuição, ou seja, a média e o desvio padrão.

Essa distribuição é uma distribuição binomial. E podemos encontrar a média e o desvio padrão de uma distribuição binomial pelas fórmulas teóricas, na qual Prob é a probabilidade de um resultado cara se a moeda é honesta.

Teoricamente a média dessa distribuição é \(n \times Prob = 1.000 \times 0.5 = 500\) (como já imaginávamos)

A variância desssa distribuição é:

\(n \times Prob \times (1- Prob) = 1000 \times 0.5 \times (1-0.5)\)

variância = 1000 x 0.5 x 0.5 = 250 desvio padrão = sqrt (250) = 15.81139

Podemos também calcular a média e o desvio padrão dessa distribuição amostral usando a distribuição que criamos com o R usando a função replicate na seção anterior. Lembre-se que colocamos esses dados na variável n.caras.

Vamos calcular a média e o desvio padrão da distribuição que criamos e colocar o resultado nas variáveis m.caras e sd.caras e então comparar o resultado que obtivemos com as fórmulas teóricas com o resultado que obtivemos com essa distribuição que criamos:

# calculando a média 
m.caras   <- mean(n.caras)
m.caras
## [1] 500.1612
# caculando o desvio padrão
sd.caras  <- sd(n.caras)
sd.caras
## [1] 15.95115

Verifique que os valores observados na distribuição construida no R foram muito similares aos valores teóricos esperados de uma distribuição binomial.

Usando esses valores teóricos da média e do desvio padrão, valos calcular o percentual de vezes nos quais a soma das caras foi maior que 545. Mas antes de fazer os cálculos, vamos verificar visualmente o que desejamos:

# parametros da curva teórica
media.teorica <- 500      # média teórica
dp.teorico    <- 15.81139 # desvio padrão teórico já calculado acima

# Gerando os pontos da curva normal na região a ser hachurada
cord.x  <- c(545,    seq(545,560,0.01),    560)
cord.y  <- c(0,dnorm(seq(545,560,0.01), mean = 500, sd = dp.teorico),0)

# plotando a curva de distribuição teórica - nossa distribuição nula
curve(dnorm(x,mean=media.teorica,sd = dp.teorico),
      xlim=c(440,560),
      main='Distribução Normal Teórica - Distribuição Nula') 

# Colorindo o poligono da região do valor de p
polygon(cord.x,cord.y,col='red')

Podemos fazer um ZOOM na área de interesse, mostrando a região acima de 545, hachurada em vermelho. Essa área representa a probabilidade de encontrarmos 545 caras, levando em conta que a moeda não é viciada (essa é nossa hipótese nula). Essa área representa nosso valor de p nesse teste de hipóteses.

# plotando a curva de distribuição teórica - nossa distribuição nula
curve(dnorm(x,mean=media.teorica,sd = dp.teorico),
      xlim=c(540,560),
      main='Distribução Normal Teórica - Distribuição Nula') 

# Colorindo o poligono da região do valor de p
polygon(cord.x,cord.y,col='red')

Vamos agora aprender a calcular essa probabilidade a partir dessa a curva normal teórica.

Usando função pnorm para calcularmos a probabilidade de um evento na curva normal teórica

Para calcular o valor de p com a função pnorm() precisamos informar nos argumentos da função o ponto de interesse (nesse caso 545), a média da curva (nesse caso media.teorica = 500), e o desvio padrão da curva (dp.teorico = 15.81139 )

1-pnorm(545,mean=media.teorica, sd=dp.teorico)
## [1] 0.002213265

ou então

pnorm(545,mean=media.teorica, sd=dp.teorico, lower.tail = FALSE)
## [1] 0.002213265

Vamos fazer agora a operação inversa, vamos usar a função qnorm(), para calcular o ponto crítico, acima do qual há 5% das observações. Lembrando que a função qnorm() também avalia a área à esquerda. Precisamos informar a àrea à esquerda do ponto.

qnorm(1-0.05, mean=media.teorica, sd=dp.teorico)
## [1] 526.0074

Existe também outra maneira de fazer esse cálculo, podemos informar que a área inserida na função é a área à direita do ponto desejado.

qnorm(0.05, mean=media.teorica, sd=dp.teorico, lower.tail = FALSE)
## [1] 526.0074

Esse resultado significa que o ponto acima marca o limite no eixo x, de forma que a área abaixo da curva à direita desse ponto é de 0.05.

Ou seja, significa que a probailidade de resultados iguais ou maiores que esse tem uma chance de ocorrer de 5% (0.05).

Discussão final

Espera-se que ao final dessa simulação o aluno tenha experimentado e compreendido as bases de uma simulação simples no R, compreendido o que é a a função de uma distribuição amostral, e os passos para o cálculo do valor de p para fins de decisão sobre a rejeição ou não da hipótese nula.

Extras - Gráficos com ggplot2:

dando um ZOOM na região de interesse. O que queremos é calcular a área abaixo da curva, à direita de 545. Essa área representa a probabilidade de encontrarm os valores iguais ou maiores que 545.