Probabilidade
Prof. Lupercio F. Bessegato
11/11/2020
Probabilidade de eventos:
Problemas de probabilidade podem ser difíceis de serem resolvidos. Uma alternativa válida é recorrer aos recursos computacionais pra auxiliar com os cálculos tediosos ou mesmo de aproximar uma solução por meio de simulações. O R apresenta funções que são úteis para simular sorteios de elementos de um conjunto finito qualquer, ou de números de acordo com padrões de distribuições de probabilidade de variáveis aleatórias. Vamos focar o primeiro caso, ou seja, sortear uma quantidade qualquer de elementos de um conjunto finito qualquer. O sorteio pode se dar com ou sem reposição e a probabilidade de se escolher um elemento desse conjunto pode ser a mesma ou não. A função sample() é muito útil para essa finalidade. Ela gera uma amostra de tamanho especificado de conjunto de dados ou de elementos, com ou sem reposição e é usada para obter amostras de vetores numéricos, de caracteres ou de data frames. Seus argumentos são:
x: vetor com os elementos a serem sorteados (ou um inteiro positivo).
size: quantidade de itens a serem sorteados.
replace: valor lógico indicando se a amostra deve ser gerada com reposição (default: FALSE).
prob: vetor com o peso para obter os elementos de x.
Entenda seus argumentos com os exemplos abaixo:
# sample() com reposição
sample(1:20, 10, replace = TRUE)## [1] 16 6 9 2 2 18 8 10 18 15sample(LETTERS, 5, replace = TRUE)## [1] "U" "G" "B" "R" "W" # sample sem reposição
sample(1:20, 10, replace = FALSE)## [1] 12 16 4 3 11 18 19 17 9 7sample(LETTERS, 5, replace = FALSE)## [1] "N" "Y" "T" "Z" "K"pode então gerar permutações, por exemplo gerar 5 permutações da palavra AMOR:
amor <- c("A", "M", "O", "R")
for(i in 1:5) print(sample(amor, 4, replace = FALSE))## [1] "M" "O" "A" "R"
## [1] "O" "M" "R" "A"
## [1] "R" "M" "A" "O"
## [1] "M" "A" "O" "R"
## [1] "R" "A" "M" "O"Se quisermos formar alguns anagramas temos de juntar as letras que são sorteadas sem reposição. Utilizamos a função paste() com o argumento collapse="". Verifique agora os resultados:
# formando palavras
for(i in 1:5) print(paste(sample(amor, 4, replace = FALSE),
collapse=""))## [1] "AROM"
## [1] "RAMO"
## [1] "MORA"
## [1] "MORA"
## [1] "ORMA"Uma outra facilidade na formação de anagramas aleatórios é utilizar a função strsplit() para transformar qualquer string em uma lista de caracteres com cada uma das letras que formam a palavra "amor":
strsplit(amor, "")## [[1]]
## [1] "A"
##
## [[2]]
## [1] "M"
##
## [[3]]
## [1] "O"
##
## [[4]]
## [1] "R"ainda, não conseguimos transformar uma string em um vetor de caracteres. O comando unlist() oferece essa funcionalidade. Confira o resultado:
# transformando string em vetor com cada um de seus caracteres
unlist(strsplit(amor, ""))## [1] "A" "M" "O" "R"Podemos agora criar uma função para formar anagrama aleatório de uma palavra sem letras repetidas ou espaços:
# função para criar anagramas da palavra x
anagrama <- function(palavra){
# transformando palavra em vetor de caracteres
letras <- unlist(strsplit(palavra, ""))
# quantidade de letras
qte <- length(letras)
# formção do anagrama
anag <- paste(sample(letras, qte, replace = FALSE), collapse="")
print(anag)
}Verifique o funcionamento da função, por exemplo, com a palavra “PESQUISA”:
anagrama("PESQUISA")## [1] "SEPQUAIS"Claro que a cada repetição da função anagrama() a palavra será formada aleatoriamente.
Podemos também formar permutações envolvendo números, como por exemplos de permutação do saudoso ano de 2019:
# permutação de 2019
numero <- sample(c(2, 0, 1, 9), 4, replace = FALSE)
as.numeric(paste(numero, collapse = ""))## [1] 1902Foi utilizada a função as.numeric(), pois o coamando paste(vetor, collapse = "") concatena os elementos de vetor em uma unica string.
Outro argumento da função sample() é prob, o qual indica a probabilidade de cada item que pode ser sorteado. Suponha uma urna com 10 bolas, com três diferentes cores: 5 delas são "red"; 2, "blue4" e 3,"chartreuse". Podemos então usar o argumento prob = c(0.5, 0.2, 0.3)
# urna com 10 bola
bolas <- c("red", "blue4", "chartreuse")
# simulação da retirada de 3 das 10 bolas
sample(bolas, 3, prob = c(0.50, 0.20, 0.30), replace = TRUE)## [1] "red" "red" "chartreuse"Veja o gráfico abaixo que simula a retirada de 10 bolas dessa urna, usando palete rainbow:
# sorteio da posição da cor
retirada <- sample(1:3, 10, prob = c(0.50, 0.20, 0.30), replace = TRUE)
# centros dos círculos
circ.x <- rep(c(2, 4, 6, 8, 10), each = 2)
circ.y <- rep(c(4, 8), times = 5)
# raios
uns <- rep(1, 10)
# urna com as bolas
symbols(x = circ.x, y = circ.y, circles = 0.5 * uns, inches = FALSE,
xlim = c(1, 11), ylim = c(1, 10),
xlab = "", ylab = "", main = "Bolas sorteadas",
xaxt = "n", yaxt = "n",
bg = bolas[retirada], fg = "gray30")
# numeração das bolas sorteadas
text(x = circ.x, y = circ.y, cex = 0.8, font = 2)
Outro uso importante da função sample() na amostragem sem reposição de observações de um conjunto de dados. Vamos amostrar 5 observações do conjunto de dados mtcars (veja sua documentação). Na realidade escolheremos ao acaso 5 das 32 linhas deste conjunto de dados.
# sorteio sem reposição de data frame
dim(mtcars)## [1] 32 11 # linhas sorteadas (nenhuma se repete, mesma probabilidade de escolha)
linhas <- sample(nrow(mtcars), 5, replace = FALSE)
linhas## [1] 15 28 9 24 3mtcars[linhas, ]## mpg cyl disp hp drat wt qsec vs am gear carb
## Cadillac Fleetwood 10.4 8 472.0 205 2.93 5.250 17.98 0 0 3 4
## Lotus Europa 30.4 4 95.1 113 3.77 1.513 16.90 1 1 5 2
## Merc 230 22.8 4 140.8 95 3.92 3.150 22.90 1 0 4 2
## Camaro Z28 13.3 8 350.0 245 3.73 3.840 15.41 0 0 3 4
## Datsun 710 22.8 4 108.0 93 3.85 2.320 18.61 1 1 4 1Quando executamos o código acima:
A função `sample()` escolhe ao acaso 5 linhas para indexar o objeto.
Os 5 índices são usados como entrada para buscar essas 5 linhas.A função sample() pode ser usada para amostrar ao acaso com reposição observações de um conjunto. Há procedimentos importantes em Estatística Computacional que usam esse recurso (por exemplo: bootstrap). Vamos sortear ao caso, com reposição 15 das 150 linhas do conjunto de dados iris.
# sorteio com reposição de data frame
indices <- sample(1:nrow(iris), 15, replace = TRUE)
indices## [1] 133 130 6 97 37 69 23 117 8 35 139 17 57 91 150head(iris[indices, ])## Sepal.Length Sepal.Width Petal.Length Petal.Width Species
## 133 6.4 2.8 5.6 2.2 virginica
## 130 7.2 3.0 5.8 1.6 virginica
## 6 5.4 3.9 1.7 0.4 setosa
## 97 5.7 2.9 4.2 1.3 versicolor
## 37 5.5 3.5 1.3 0.2 setosa
## 69 6.2 2.2 4.5 1.5 versicolorAgora é possível a repetição de linhas. Compare a a aplicação da função prop.table() tanto ao conjunto completo quanto à essa amostra escolhida com reposição:
# proporção de cada espécie - cjto completo
prop.table(table(iris[, 5]))##
## setosa versicolor virginica
## 0.3333333 0.3333333 0.3333333 # proporção de cada espécie - amostra com reposição
prop.table(table(iris[indices, 5]))##
## setosa versicolor virginica
## 0.4000000 0.2666667 0.3333333O que aconteceria com essas porcentagens se repetíssemos esse procedimento um grande número de vezes? Para isso, podemos usar a função replicate() que repete uma determinada função por uma quantidade especificada de vezes
Por exemplo, para repetir as 5 simulações do lançamento de 1 moeda, 8 vezes (considere como resultados: 0, representando cara e 1, coroa);
# 5 simulações de 8 lançamentos
replicate(5, sample(0:1, size = 8, replace = TRUE))## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0 0 0 0
## [2,] 0 1 0 1 0
## [3,] 1 0 0 1 1
## [4,] 1 0 1 1 0
## [5,] 0 0 0 0 0
## [6,] 0 1 1 0 1
## [7,] 1 1 1 0 0
## [8,] 1 0 1 0 0O resultado é uma matriz cujo número de lihas é a quantidade de simulações (repetições dos 8 lançamentos) e a quantidade de colunas é o tamanho da função que foi repetida (no exemplo um vetor com 8 valores de 0’s ou 1’s). vamos montar agora uma simulação da proporção de 1’s em 8 lançamentos de moedas (a média de vetor de 0’s e 1’s):
# proporção de 1's em 10 lançamentos
replicate(5, mean(sample(0:1, size = 8, replace = TRUE)))## [1] 0.875 0.750 0.625 0.500 0.625O resultado agora é um vetor de tamanho 5, pois o resultado da função é um único valor a média. Poderíamos, visualizar o histograma da proporção de 1’s em 1.000 simulações do lançamento de uma moeda 25 vezes.
# histograma com 1.000 simulações da proporção de 1's em 25 lançamentos
hist(replicate(1000, mean(sample(0:1, size = 25, replace = TRUE))),
freq = F, main = "Proporção de 1's",
xlab = expression(hat(p)), ylab = "Densidade")
mtext("n = 25", line = 0.25, side = 3, cex = 0.85, font = 2)
Problema do aniversário:
Um problema interessante em probabilidade é o paradoxo do aniversário. Ele diz respeito à probabilidade de que, em um conjunto de n pessoas escolhidas ao acaso, pelo menos duas delas façam aniversário no mesmo dia. Leia mais a respeito deste problema muito estudado em Probabilidade: Birthday problem. Para simplificar o problema faremos algumas considerações: (i) não ha alunos nascidos os anos bissextos; (ii) não há gêmeos na turma; (iii) não há variações sazonais nas datas de nascimento, ou seja, os 365 dias de aniversários são igualmente prováveis. Isso não é real! Leia essa matéria da BBC Brasil, de 18/02/2019: Brasileiros nascem mais entre março e maio, mas razão intriga cientistas. Considerem uma turma com 16 alunos e que todos os dias de um ano não-bissexto são igualmente prováveis. É relativamente fácil de calcular a probabilidade de todos os aniversários serem diferentes nessa turma:
\[\begin{align} \operatorname{P}\{\text{16 dias são distintos} \}& = \overbrace{\frac{365}{365}\frac{364}{365}\frac{363}{365} \dots \frac{350}{365}}^{\text{16 datas}} \\ & = \frac{1}{365^{16}}365 (365 -1)(365 - 2) \dots (365 - (16-1)) \\ & = \frac{1}{365^{16}} \frac{365!}{(365-16)!} \\ \end{align}\]
A probabilidade de pelo menos uma coincidência de dias de aniversário é agora facilmente calculada:
\[\begin{align} \operatorname{P}\{\text{pelo menos uma coincidência} \}& = 1 - \frac{1}{365^{16}} \frac{365!}{(365-16)!} \end{align}\]
que pode ser induzida para uma turma com n alunos como sendo
\[\begin{align} \tag{*} \operatorname{P}\{\text{pelo menos uma coincidência em n alunos} \} & = 1 - \frac{1}{365^{n}} \frac{365!}{(365-n)!} \end{align}\]
Agora é só calcular para nosso caso, em que \(n=16\). Há, entretanto alguns problemas séries: os fatoriais e as potências são bastante elevados e pode haver problemas de exceder o máximo valor possível no R. Calcule, por exemplo, \(365!\):
factorial(365)## [1] InfE, curiosidade, qual é o maior fatorial que consigo calcular? e potência de 365? Veja abaixo:
factorial(170:171)## [1] 7.257416e+306 Inf365^(120:121)## [1] 2.986371e+307 InfTemos algumas alternativa com a fórmula \(\tag{*}\) para \(n = 16\):
1 - podemos usar \prod() para calcular a razão entre os fatoriais
# usando produtorio
1- (1/365^16)*prod(365:350)## [1] 0.283604que é a resposta desejada.
2- O coeficiente binomial de 365 por n é: \[\begin{align} \binom{365}{n} = \frac{365!}{n!(365 - n)!} \end{align}\]
logo, a equação \(\tag{*}\) pode ser reescrita como:
\[\begin{align} \operatorname{P}\{\text{pelo menos uma coincidência em n alunos} \} & = 1 - \frac{n!}{365^{n}} \binom{365}{n} \end{align}\]
O comando choose(365, 16) calcula \(\binom{365}{16}\), logo, a probabilidade desejada pode ser calculada por:
# usando choose()
1-(factorial(16)/365^16)*choose(365, 16)## [1] 0.2836043- A função lfactorial() calcula o log do fatorial desejado. Pode-se então calcular o log da razão da equação \(\tag{*}\): \[\begin{align}
\log\left( \frac{1}{365^{n}} \frac{365!}{(365-n)!} \right) =
\log(365!) - \log((365-n)!) - n \log{365}
\end{align}\]
o resultado desejado pode ser obtido por meio do exponencial desse logaritmo:
# usando lfactorial
1 - exp((lfactorial(365) - lfactorial(365-16)) - 16*log(365))## [1] 0.2836044 - Similarmente a função lgamma() calcula o logaritmo da função gama e sabemos que, se \(x \in \mathbb{N}\), então \(\Gamma(x) = (x-1)!\), assim a razão da equação \(\tag{*}\) pode ser calculada por: \[\begin{align}
\log\left( \frac{1}{365^{n}} \frac{365!}{(365-n)!} \right) =
\log(\Gamma(365 +1)) - \log(\Gamma(365-n + 1)) - n \log{365}
\end{align}\]
O resultado desejado também pode ser obtido pela função inversa desse logaritmo:
# usando lgamma
1 - exp((lgamma(365 + 1) - lgamma(365-16 + 1)) - 16*log(365))## [1] 0.2836045- Há uma função do R Base que efetua o mesmo cálculo, pbirthday():
# função built-in do R
pbirthday(16)## [1] 0.283604Agora, surge uma pergunta: qual o comportamento dessa probabilidade para outros valores de n? Vamos construir a função aniversario() para calcular essa probabilidades de uma maneira mais fácil (utilizaremos a alternativa 2 acima):
# Função para cálculo da probabilidade de dois ou mais estudantes
# ter o mesmo dia de aniversário
aniversario <- function(N){
resultado <- 1-(factorial(N)/365^N)*choose(365,N)
resultado
}O resultado para turma de 16 alunos é facilmente obtido:
aniversario(16)## [1] 0.283604Podemos agora obter um gráfico dessa probabilidade de pelo menos uma coincidência de dias de aniversário e o tamanho da turma:
# gráfico
plot(n <- 0:100, aniversario(n), type = "l", xaxs = "i", yaxs = "i",
xlab = "n", ylab = "P(pelo menos um par)",
main = "Problema do Aniversário",
ylim = c(0, 1))
mtext(expression("Função de Probabilidade"), line = 0.25,
cex = 1.1)
# n = 16
arrows(16, 0, 16, aniversario(16), length = 0.1, angle = 20,
col = "blue")
arrows(16, aniversario(16), 0, aniversario(16), length = 0.1,
angle = 20, col = "blue")
# anotação da probabilidade
text(0, aniversario(16) + 0.04, round(aniversario(16), 3),
cex = 0.7, adj = c(0, 0), offset = 1.5, font = 2, col = "blue")
Inspecionando o gráfico obtemos alguns valores ao longo dos diverssos tamanhos de turmas:
| n | Probabilidade |
|---|---|
| 10 | 0.1169482 |
| 16 | 0.2836040 |
| 20 | 0.4114384 |
| 30 | 0.7063162 |
| 40 | 0.8912318 |
| 50 | 0.9703736 |
| 100 | 0.9999997 |
e o paradoxo surge, pois, com turmas relativamente pequenas é bastante alta a probabilidade de haver pelo menos um par de estudantes que comapratilham a mesma data de aniversário. Com 100 alunos a probabilidade é bastante altísssima, para 300 alunos a probabilidade é de \(1 - 7\times 10^{-73}\).
Uma outra abaordagem é aproximar o resultado por meio de simulação. Há problemas de probabilidade que são muito difíceis de serem resolvidos analiticamente, sendo bastante conveniente aproximação empírica.
# ------------ Simulação - 1 turma de 16
# dias de aniversario de uma turma de 16 alunos (1 a 365)
aniversarios <- sample(x = 1:365, size = 16, replace = TRUE)
# dias diferentes no conjunto de datas sorteadas
unique(aniversarios)## [1] 89 334 194 285 57 41 157 55 336 148 258 361 234 99 284 18 # quanidade de datas diferentes
length(unique(aniversarios))## [1] 16 # se qte datas diferentes for menor que tamanho da turma
coincide <- length(unique(aniversarios)) < 16
# se 1: pelo menos uma coincidência se 0, não há datas iguais
sum(coincide)## [1] 0O código acima simula uma turma. Se sum(coincide) for igual a 1, houve um empate de datas de aniversário, se 0, não houve. Agora fica simples repetir o procedimento milhares de vezes, verificando a quantidade de turmas em que houve coincidência de datas. A proporção de turmas com coincidência de datas é uma aproximação da probabilidade desejada. Antes, para facilitar a simulação, será construída a função simula() para codificar o procedimento acima:
simula <- function(n){
datas <- sample(x = 1:365, size = n, replace = TRUE)
repete <- length(unique(datas)) < n
sum(repete)
}
# uma particular turma
simula(16)## [1] 0Agora, simulam-se 1000 turmas, claculando-se a proporção de turmas em que houve coincidência:
vezes <- 1000
prob.emp <- sum(replicate(vezes, simula(16)))/vezes
prob.emp## [1] 0.292O valor apresenta um erro de 0.008396 com relação ao verdadeiro valor da probabilidade desejado(0.283604). E se quisermos construir uma distribuição de probabilidade para as variás possibilidades, ou sejam, desde nenhuma coincidência até os 16 aniversários ocorrerem na mesma data.
A função duplicated() tem como saída um vetor lógico indicando quais valores são duplicados com relação a linhas com menore índices. Veja, por exemplo, o que acontece se simulamos em uma turma de 40 estudantes:
# amostra com dias de aniversário de turma de tamanho n
n <- 40
amostra <- sample(sample(x = 1:365, size = n, replace = TRUE))
repetidos <- duplicated(amostra)
# Datas que se repetem
unique(amostra[repetidos])## [1] 35 # quatidade de alunos com mesmo dia de aniversário
table(amostra[repetidos]) + 1##
## 35
## 2Uma alternativa com um código mais curto e mais elegante é usando também a função any(vetor), que fornece a resposta TRUE se há pelo menos um valor TRUE em um conjunto de vetores lógigos. (A função complementar é all(vetor), que fornece a resposta TRUE apenas se são TRUE todos os valores em um conjunto de vetores lógicos). acompanhe a execução dos comandos abaixo aplicados ao conjunto de datas simuladas do objeto aniversarios:
# --- simulação - 1 turma de 16
# (usando any())
# simulacao de turma anterior
aniversarios## [1] 89 334 194 285 57 41 157 55 336 148 258 361 234 99 284 18 # datas duplicadas - varia de 0 a 15
duplicated(aniversarios)## [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE FALSE FALSE # datas que coincidem
aniversarios[duplicated(aniversarios)]## integer(0) # TRUE: pelo menos uma coincidencia
# FALSE: nenhuma coincidencia
any(duplicated(aniversarios))## [1] FALSE # quantidade de datas que coincidem
sum(any(duplicated(aniversarios)))## [1] 0A função simula() pode ser reescrita com essa codificação:
# ---- Função para simular datas de uma turma com n pessoas
# Alternativa 2 - usando any()
simula <- function(n){
datas <- sample(x = 1:365, size = n, replace = TRUE)
sum(any(duplicated(datas)))
}e, com essa função simula(), obtemos mais uma estimativa para a probabilidade de pelo menos uma coincidência com os comandos abaixo:
# simulação de uma turma com 16 pessos
simula(16)## [1] 0 # --- Simulação d run turmas com n pessoas
run <- 10000
n <- 16
coincidencias <- sum(replicate(run, simula(n)))
# qte de turmas com pelo menos 1 coincidencia
coincidencias## [1] 2813 # proporção de pelo menos uma coincidência -
# Aproximação P(pelo menos um par de dias de aniversário iguais)
coincidencias/run## [1] 0.2813A estimativa apresenta um erro de -0.002304 com relação ao verdadeiro valor da probabilidade desejado(0.283604).
Problema de estimação de probabilidade de aprovação:
Este tópico está baseado no post The “probability to win” is hard to estimate….
Em geral, é difícil estimar em tempo real a probabilidade de ganhar como, por exemplo, em jogos de futebol, em eleições, etc. Mas, podemos ver isso em apostas com questões de múltipla escolha.
Considere uma avalição de múltipla escolha clássica. Após cada pergunta, imagine que você queira calcular a probabilidade de o aluno ser aprovado. Considere aqui o caso em que temos 50 questões. Os alunos são aprovados quando têm pelo menos 25 respostas corretas. Para simular a aplicação dessa avaliação, considere que os estudantes decidem a resposta de cada pergunta jogando uma moeda. São \(n\) alunos e 50 perguntas por avaliação, Considere \(n=10\) e armazene as respostas na matriz de zeros e uns \(M_{50\times 10}\):
# quantidade de alunos
n <- 10
# matrix com as respostas
M <- matrix (sample (0: 1, size = n * 50, replace = TRUE), nrow = 50)
head(M)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,] 1 0 0 0 0 1 0 0 0 0
## [2,] 0 0 1 1 0 0 0 1 0 1
## [3,] 0 0 1 1 0 0 0 0 0 0
## [4,] 1 1 0 0 1 1 1 1 1 0
## [5,] 0 0 0 1 0 1 1 0 0 0
## [6,] 0 0 0 0 1 1 0 0 0 0Seja \(X_{i, j}\) a pontuação do aluno \(j\) na questão \(i\) e \(S_{i, j}\), a pontuação acumulada, ou seja, \(S_{i, j} = X_{1, j} + \dots + X_{i, j}\). Na etapa $ i $, poderia obter algum tipo de previsão da pontuação final, usando \(\hat{T}_{i, j} = 50 \times S_{i, j}/i\). O código abaixo implementa o estimador \(\hat{T}_{i, j}\):
# matriz com soma acumulada por coluna
SM <- apply(M,2,cumsum)
# predição de acertos após a correção de cada questão
NB <- SM*50/(1:50)SM e NM são, respectivamente, as matrizes dos resultados acumulados das correções e da predição do resultado final baseada nessa soma acumulada. Veja por exemplo a predição do estimador da pontuação dos alunos (todas as colunas) \(\hat{T}_{25, j}\), ao término da correção de metade da avaliação (NB[25,]) e o resultado final (SM[50, ])
# predição após a correção de metade da prova
NB[25, ]## [1] 18 22 22 24 20 26 22 18 22 18 # ultima linha - qte de acertos de cada um dos 10 alunos
SM[50, ]## [1] 21 24 23 22 25 26 21 25 26 21Verifique o gráfico para os alunos ("lightblue") e para o 3º aluno ("red")
# predição do aluno 1 durante a correção da avaliação
plot(NB[, 1], type = "n", ylim = c(0, 50))
# linha limite entre aprovação e reprovação
abline(h = 25, col = "blue")
# curva das predições dos outros alunos
for(i in 1:n) lines(NB[, i], type ="s", col = "light blue")
# linha de predição do aluno 3
lines(NB[, 3], type = "s", col = "red")
Mas estamos simplesmente predizendo, em cada etapa, a pontuação final. Essa não é estimação da probabilidade de aprovação! Há situações a serem consideradas: (i) Se após \(i\) questões, \(i \geq 25\), o aluno \(j\) possuir 25 respostas corretas, é trivial estabelecer que a probabilidade de aprovação é 1, ou seja, se $S_{i, j} geq 25 $; (ii) outro caso simples para determinação da probabilidade de aprovação ocorre quando, após as perguntas \(i\), o número de pontos que ainda é possivel ele obter até o final for insuficiente para alcançar o total de 25 acertos necessários para a aprovação, ou seja, se \(S_{i, j} + (50 - i + 1) < 25\) a probabilidade de aprovação do aluno \(j\) é 0; (iii) nas demais situações, a probabilidade de aprovação do aluno \(j\) é simples de ser calculada: é a probabilidade de obter pelo menos \(25 - S_{i, j}\) respostas corretas, nas \(50 - i\) questões faltantes, considerando como probabilidade de sucesso \(S_{i, j}/i\) (considerada aqui como um avaliador de seu conhecimento). Essa é a função de sobrevivência (\(1 - F_X(x)\)) de uma distribuição binomial. O código é simples de ser implementado;
# ---------- Probabilidade de aprovação do aluno j,
# dado o desempenho em i questões
# matriz das probabilidades de aprovação
PB=NB*NA
# regra para cálculo das probabilidades
for(i in 1:50){
for(j in 1:n){
if(SM[i,j] >= 25) PB[i,j] = 1
if(SM[i,j] + (50 - i + 1) < 25) { PB[i,j] = 0
}else{
PB[i,j] = 1- pbinom(25 - SM[i, j], size = (50- i),
prob = SM[i, j]/i)
}
}
}
tail(PB)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## [45,] 0 0.1739197 0 0 0.02213267 0.1739197 0 0.02213267 0.5208265
## [46,] 0 0.3125000 0 0 0.05231900 0.3125000 0 0.05231900 0.7193799
## [47,] 0 0.1171898 0 0 0.11718983 0.5159550 0 0.11718983 0.5159550
## [48,] 0 0.0000000 0 0 0.25000000 0.2500000 0 0.25000000 0.7703993
## [49,] 0 0.0000000 0 0 0.51020408 0.5102041 0 0.51020408 0.5102041
## [50,] 0 0.0000000 0 0 0.00000000 1.0000000 0 0.00000000 1.0000000
## [,10]
## [45,] 0
## [46,] 0
## [47,] 0
## [48,] 0
## [49,] 0
## [50,] 0Podemos plotar essa matriz de probabilidades de aprovação, estimadas após a correção de cada questão (correção ordenada):
# --------- Plot da matriz de probabilidades de aprovação
# Estimação após cada correção
plot(PB[, 1], type = "n", ylim = c(0,1))
for(i in 2:n) lines(PB[, i], type = "s", col = "lightblue")
lines(PB[, 3], type = "s", col = "red")
As estimativas de probabilidade de aprovação são menos voláteis que as estimativas de acerto. Perceba que as probabilidade de aprovação convergem ou para 1 (aprovação) ou para 0 (reprovação). Claro que o comportamento será ligeiramente diferente se os alunos não usarem uma moeda para responder cada questão, mas se consideramos, por exemplo, que metade dos alunos dominam o conteúdo (probabilidade de acerto da questão igual a 2/3) ou que não tiveram tempo de estudar o conteúdo (probabilidade de acerto igual a 1/3). A codificação será a mesma, exceto para a construção da matriz de respostas M:
# ------- Probabilidades de acerto diferentes ----
#
# quantidade de alunos
n <- 10
# matrix com as respostas
M <- matrix (sample (0:1 , size = n * 50, prob = c(1/3, 2/3),
replace = TRUE), nrow = 50)
# troca de sucessos por fracassos (probabilidades complementares)
M[, 6:10] <- 1 - M[, 6:10]
# quantidade de acertos de cada um dos 10 alunos
colSums(M)## [1] 35 33 36 33 36 21 21 11 8 18 # --------- Predição Móvel ------------
# matriz com soma acumulada por coluna
SM <- apply(M, 2, cumsum)
# ultima linha - qte de acertos de cada um dos 10 alunos
SM[50, ]## [1] 35 33 36 33 36 21 21 11 8 18 # predição de acertos após a correção de cada questão
NB <- SM*50/(1:50)
# predição após a correção de metade da prova
# ---------- Probabilidade de aprovação do aluno j,
# dado o desempenho em i questões
# matriz das probabilidades de aprovação
PB=NB*NA
# regra para cálculo das probabilidades
for(i in 1:50){
for(j in 1:n){
if(SM[i,j] >= 25) PB[i,j] = 1
if(SM[i,j] + (50 - i + 1) < 25) { PB[i,j] = 0
}else{
# if((SM[i,j] < 25) & (SM[i, j] + (50 - i + 1) >= 25))
PB[i,j] = 1- pbinom(25 - SM[i, j], size = (50- i),
prob = SM[i, j]/i)
}
}
} # ----------- Gráfico da predição móvel ----
# predição do aluno 1 durante a correção da avaliação
plot(NB[, 1], type = "n", ylim = c(0, 50))
# linha limite entre aprovação e reprovação
abline(h = 25, col = "blue")
# curva das predições dos outros alunos
for(i in 1:n) lines(NB[, i], type ="s", col = "light blue")
# linha de predição do aluno 3
lines(NB[, 3], type = "s", col = "red")
A estimação das probabilidades apontam que, na maioria dos casos, não é necessário aguardar até o final da correção para saber quem foi aprovado ou não.
# --------- Plot da matriz de probabilidades de aprovação
# Estimação após cada correção
plot(PB[, 1], type = "n", ylim = c(0,1))
for(i in 2:n) lines(PB[, i], type = "s", col = "lightblue")
lines(PB[, 3], type = "s", col = "red")
Talvez as estimações apresentem menor variabilidade caso se adote uma abordagem Bayesiana para o problema.
Distribuições de probabilidade no R:
O R Básico tem implementada a maioria das distribuições de probabilidade. Iremos mostrar como utilizar as funções relacionadas com essas distribuições, com foco em sua visualização. Isso auxiliará no entendimento de como as diferentes funções do R estão efetuando os cálculos.
Cada distribuição de probabilidade está associada com quatro funções que seguem a seguinte convenção para seus prefixos:
d: função de densidade de probabilidade, se contínua ou função de probabilidade (mass function), se discreta da distribuição de probabilidade.
p: função de distribuição de probabilidade acumulada.
q: inversa da função de distribuição acumulada, ou função quantil da distribuição de probabilidade.
r: gerador de números aleatórios da distribuição de probabilidade.
Para efetuar o cálculo desejado, cada função usa um argumento, x, q, p e n, respectivamente, seguido dos parâmetros da distribuição que definem a particular distribuição de probabilidade.
| Distribuição | Núcleo | Parâmetros | Default |
|---|---|---|---|
| Beta | beta |
shape1,shape2 |
|
| Binomial | binom |
size, prob |
|
| Binomial negativa | nbinom |
size,prob,mu |
|
| Cauchy | cauchy |
location,scale |
0, 1 |
| Exponencial | exp |
1/mean |
1 |
| F | f |
df1, df2 |
|
| Gama | gamma |
shape,1/scale |
NA, 1 |
| Geométrica | geom |
prob |
|
| Hipergeométrica | hyper |
m, n, k |
|
| Log-normal | lnorm |
mean, sd |
0, 1 |
| Logística | logis |
location,scale |
0, 1 |
| Normal | norm |
mean, sd |
0, 1 |
| Poisson | pois |
lambda |
|
| Qui-quadrado | chisq |
df |
|
| t de Student | t |
df |
|
| Uniforme | unif |
min, max |
0, 1 |
| Weibull | weibull |
shape |
A seguir, apresentamos alguns gráficos e exemplos de uso da função de densidade e da função quantil.
Distribuição normal:
A distribuição normal \(\mathcal{Normal}(\mu, \sigma)\), onde \(\mu\) é a média e \(\sigma\), o desvio-padrão é representada no R pelas funções dnorm, pnorm e qnorm, A função de densidade de probabilidade dnorm() e a função de distribuição acumulada, pnorm() são definidas em toda a reta. O parâmetro \(\mu\) é um parâmetro de locação pois controla a posição de seu ponto de simetria, por outro lado, \(\sigma\) é seu parâmetro de escala, pois controla a dispersão da massa de probabilidade em torno de seu centro (\(\mu\)). A forma da distribuição normal é a mesma, não importam os valores dos parâmetros. Perceba como o valor do parâmtero de locação (\(\mu\)) influencia o gráfico da função de densidade de probabilidade da distribuição normal:
# Parâmetros das distribuições
medias <- c(5, 10, 15)
desvio <- 2
intervalo.x <- 3.5*2*c(-1, 1) + range(medias)
intervalo.y <- c(0, dnorm(0, sd = desvio) + 0.01)
cores <- c("black", "red", "blue")
# gráfico primário
plot(intervalo.x, type = "n",
xlim = intervalo.x, ylim = intervalo.y,
ylab = "Densidade", xlab = "X",
main = "Função de Densidade de Probabilidade")
# gráficos secundários e anotações das distribuições
for (i in 1:3) {
# gráfico da normal de medias[i]
curve(dnorm(x, mean = medias[i], sd = desvio),
from = medias[i] - 3.5*desvio, to = medias[i] + 3.5*desvio,
col = cores[i], lty = i,
add = TRUE)
# anotação do valor da média
text(medias[i], 0.10, bquote(mu==.(medias[i])), cex = 0.8,
col = cores[i], lty = i)
}
# conclusão do título do gráfico (2a. linha)
mtext(expression("Distribuição normal - ("~ mu ~"," ~sigma==2~")"),
cex = 0.9)
Agora, perceba a influência do parâmetro de escala (\(\sigma\)) na dispersão da função de densidade de probabilidade da distribuição normal:
# Parâmetros das distribuições
media <- 10
desvios <- c(1, 2, 4)
intervalo.x <- 3.5*max(desvios)*c(-1, 1) + media
intervalo.y <- c(0, dnorm(0, sd = min(desvios)) + 0.01)
# gráfico primário
plot(intervalo.x, type = "n",
xlim = intervalo.x, ylim = intervalo.y,
ylab = "Densidade", xlab = "X",
main = "Função de Densidade de Probabilidade")
# gráficos secundários e anotações das distribuições
i <- 1
for (i in 1:3) {
# gráfico da normal de medias[i]
curve(dnorm(x, mean = media, sd = desvios[i]),
from = media - 3.5*desvios[i], to = media + 3.5*desvios[i],
col = cores[i], lty = i,
add = TRUE)
text(media + 1.2 * desvios[i], dnorm(0, sd = desvios[i]), bquote(sigma ==.(desvios[i])), cex = 0.8,
col = cores[i], lty = i)
}
# conclusão do título do gráfico (2a. linha)
mtext(expression("Distribuição normal - ("~ mu == 10~"," ~sigma ~")"),
cex = 0.9)
Vamos verificar agora a distribuição normal padrão \(Z\), ou seja, com $= 0 $ e \(\sigma = 1\). A densidade é muito pequena fora do intervalo \((-3,5; 3,5)\), portanto o gráfico de sua função de densidade de probabilidade será restrito a esse intervalo.
# parâmetros da distribuição
mi <- 0
dp <- 1
# gráfico da densidade
curve(dnorm(x, mean = mi, sd = dp), from = -3.5, to = 3.5,
ylab = "Densidade",
main = "Função de Densidade de Probabilidade")
mtext(expression("Distribuição normal - ("~ mu==0 ~"," ~sigma==1~")"),
cex = 0.9)
lines(c(-3.5, 3.5), c(0,0))
# máxima densidade
(dens.max <- dnorm(0, mean = 0, sd =1)) ## [1] 0.3989423points(0, dens.max, pch = 21, bg = "red")
# valor da máxima densidade
text(-0.5, dens.max, paste("(0, ",round(dens.max, 3),")"), cex = 0.75,
pos = 2, adj = 1)
# seta indicativa
arrows(-0.6, dens.max, 0, dens.max, length = 0.10, angle = 15)
# probabilidade abaixo de x1
z1 <- -1.5
pnorm(z1, mean = mi, sd = dp)## [1] 0.0668072 # área da probabilidade
base <- seq(-3.5, z1, by = 0.01)
polygon(x = c(base, -1.5, -3.5),
y = c(dnorm(base, mean = 0, sd =1), 0, 0),
col = "gray86", border = "gray46")
# anotação probabilidade no gráfico
text(-2.5, 0.1, round(pnorm(-1.5), 3), cex = 0.85)
# seta indicativa
arrows(-2.5, 0.1 - 0.015, -2, dnorm(-2), length = 0.10, angle = 15)
A função de distribuição acumulada é definida como: \[ F_Z(x) = \operatorname{P}\{Z \leq x \}. \] logo, para qualquer variável aleatória, ela é crescente, com domínio em toda a reta e imagem no intervalo \([0, 1]\). No caso da normal (e das variáveis contínuas) a função de distribuição acumulada é contínua. Os valores das tabela da normal padrão (bastante usada em sala de aula) são o valores obtidos na função de distribuição acumulada, que por sua vez estão relacionados com as áreas correspondentes da função de densidade de probabilidade. Verifique no gráfico da função de distribuição acumulada da distribuição normal padrão abaixo:
# parâmetros da distribuição
mi <- 0
dp <- 1
# gráfico da função de distribuição acumulada
curve(pnorm(x, mean = mi, sd = dp), from = -3.5, to = 3.5,
xaxs = "i", yaxs = "i",
ylab = "Densidade",
main = "Função de Distribuição Acumulada")
mtext(expression("Distribuição normal - ("~ mu==0 ~","~ sigma==1~")"),
cex = 0.9)
# Valor para x = -1.5 e x = 1.5
xis <- c(-1.5, 1.5)
(F.xis <- pnorm(xis, mean = mi, sd = dp))## [1] 0.0668072 0.9331928 # setas indicativas
arrows(xis[1:2], 0, xis[1:2], F.xis[1:2], length = 0.10, angle = 20,
col = "red")
arrows(xis[1:2], F.xis[1:2], -3.5, F.xis[1:2], length = 0.10,
angle = 20, col = "red")
# valor no eixo y
mtext(round(F.xis, 4), at = F.xis, line = 0.25, side = 2, las = 1,
cex = 0.7, font = 2, col = "red")
Veja no gráfico o valor correspondente a \(x_1 = -1,5\) e \(x_2 = 1.5\), ou seja \(F_Z(x_1) = \operatorname{P}\{Z \leq x_1\}= 0.0668072\) e \(F_Z(x_2) = \operatorname{P}\{Z \leq x_1\} = 0.9331928\). Com a função de distribuição acumulada podemos calcular a probabilidade em qualquer intervalo da variável aleatória \(Z ~ \mathcal{Normal(0,1)}\), como por exemplo ${-1,5 Z 1,5 } que pode ser calculado por \(F_Z(1.5) - F_Z(-1.5)\), ou seja:
pnorm(1.5, mean = 0, sd =1) - pnorm(-1.5, mean = 0, sd =1)## [1] 0.8663856Poderíamos tomar o problema inverso, ou seja dada uma probabilidade \(p\), queremos determinar o valor \(x\), tal que \(F(x) = \operatorname{P}\{Z \leq x \} = p\). Como a função de distribuição acumulada da normal padrão é estritamente crescente, mais formalmente, queremos determinar \(F_Z^{-1}(p) = x\). Vamos tentar determinar os valores \(x_1\) e \(x_2\), tal que \(\operatorname{P} \{ x_1 \leq Z \leq x_2\} = 0,95\) de maneira simétrica, ou seja, \(\operatorname{P} \{ Z \leq x_1 \} = \operatorname{P} \{ Z \geq x_2 \} = 0,025\). É importante perceber que desejamos determinar \(x_1 = F_Z^{-1}(0,025)\) e \(x_2 = F_Z^{-1}(0,975)\).
curve(pnorm(x, mean = mi, sd = dp), from = -3.5, to = 3.5,
xaxs = "i", yaxs = "i",
ylab = "Densidade",
main = "Função de Distribuição Acumulada")
mtext(expression("Distribuição normal - ("~ mu==0 ~","~ sigma==1~")"),
cex = 0.9)
# Valor para p = 0.025 e p = 0.975
p <- c(0.025, 0.975)
(quantis <- qnorm(p, mean = mi, sd = dp))## [1] -1.959964 1.959964# setas indicativas
arrows(-3.5, p[1:2], quantis[1:2], p[1:2], length = 0.10, angle = 20,
col = "blue")
arrows(quantis[1:2], p[1:2], quantis[1:2], 0, length = 0.10,
angle = 20, col = "blue")
# valores no eixo y
mtext(round(p, 3), at = p, line = 0.25, side = 2, las = 1,
cex = 0.7, font = 2, col = "blue")
# valores no eixo x
mtext(c(expression(x[1]), expression(x[2])), at = quantis,
side = 1, cex = 0.9, col = "blue", adj = c(0, 1))
Veja no gráfico o valor correspondente a \(p_1 = 0,025\) e \(p_2 = 0.975\), ou seja \(F^{-1}_Z(0,025) = -1.959964 \Rightarrow \operatorname{P}\{Z \leq -1.959964\}= 0.025\) e \(F^{-1}_Z(p_2) = 1.959964 \Rightarrow \operatorname{P}\{Z \leq 1.959964\} = 0.975\). Dessa maneira \(\operatorname{P}\{-1.96 \leq Z \leq 1.96 \}= 0,95\). A inversa da função de distribuição acumulada, \(F^{-1}_Z\), é denominada função quantil, com domínio no intervalo \([0,1]\), tendo como imagem a reta. Seu gráfico é:
curve(qnorm(x, mean = mi, sd = dp), from = 0, to = 1,
xlab = "p", ylab = "Quantil",
main = "Função Quantil")
mtext(expression("Distribuição normal - ("~ mu==0 ~","~ sigma==1~")"),
cex = 0.9)
Exemplos:
Suponha que a variável aleatória \(X\) tenha distribuição normal com média \(\mu = 6\) e desvio-padrão \(\sigma = 4\). Determine: ${ -1 X 3 } $
pnorm(3, 6, 4) - pnorm(-1, 6, 4)## [1] 0.1865682Suponha que a variável aleatória \(X\) esteja distribuída normalmente com média \(\mu = 20\) e desvio-padrão \(\sigma = 10\). Determine o intervalo com 90% de confiança em torno da média para o valor esperado de \(X\):
sup <- qnorm(0.95, 20, 10)
inf <- qnorm(0.05, 20, 10)
c(sup, inf)## [1] 36.448536 3.551464Distribuição de Poisson:
A variável aleatória discreta \(X \sim \mathcal{Poisson}(\lambda)\), definida para os números inteiros não-negativos é representada em R por: sua função de probabilidade(mass function), dpois(), sua função de probabilidade acumulada, ppois() e sua função quantil, qpois(). Seja por exemplo \(\lambda = 2.5\).
# parâmetros da distribuição
taxa <- 2.5
x.max <- qpois(0.999, lambda = taxa)
xis <- 0:x.max
dpois(2, 2.5)## [1] 0.2565156# gráfico da densidade
plot(xis, dpois(xis, taxa), type='h',
yaxs = "i", #ylim = c(0, 0.26),
xlab = "x", ylab = expression(f[X](x)),
main = "Função de Probabilidade")
points(xis, dpois(xis, taxa), pch = 20, xpd = TRUE)
mtext(expression("Distribuição de Poisson - ("~ lambda== 2.5 ~")"),
cex = 0.9)
A função de probabilidade é representada apenas em pontos, pois, toda a massa de probabilidade está localizada nos números inteiros não-negativos, ou seja $ f_X(x) = {X = x }$ e \(\sum_{x=0}^\infty f_X(x) =1\). No exemplo, o valor mais provável (moda) é \(X=2\), com probabilidade de ocorrência:
dpois(2, lambda = taxa)## [1] 0.2565156A distribuição de Poisson apresenta estatísticas descritivas com algumas características interessantes:
A média e a variância da variável aleatória são iguais ao paraêmetro, ou seja, \(\operatorname{E}(X)=\operatorname{Var}(X) = \lambda\).
O coeficiente de variação é \(\frac{1}{\sqrt{\lambda}}\).
Se \(\lambda\) é não-inteiro, então o valor mais provável da variável aleatória \(X\), ou seja, sua moda é igual a \(\lfloor \lambda \rfloor\), o maior inteiro menor ou igual a \(\lambda\) (em R,
floor()). Se \(\lambda\) é inteiro, então as modas são \(\lambda\) e \(\lambda -1\).A mediana (\(\nu\)) está limitada a \(\lambda -\ln 2\leq \nu <\lambda + {\frac {1}{3}}\), um intervalo relativamente curto.
Tente observar essas propriedades nos próximos gráficos e nos exercícios.
No caso da Poisson (e das variáveis discretas) a função de distribuição acumulada é descontínua, já que a massa de probabilidadese concentra em ponto do contradomínio da variável aleatória. Verifique o gráfico da função de distribuição acumulada da distribuição de Poisson com parâmetro \(\lambda = 2,5\) abaixo:
# Parâmetros do gráfico
taxa <- 2.5
x.max <- qpois(0.999, lambda = taxa)
xis <- 0:x.max
acumulada <- ppois(xis, taxa)
# função de distribuição acumulada - sem pontos
plot(xis, acumulada, type='n',
#xaxs = "i",
yaxs = "i", ylim = c(0, 1),
xlab = "x", ylab = expression(F[X](x)),
main = "Função de Distribuição Acumulada")
# 2a. linha do título do gráfico
mtext(expression("Distribuição de Poisson - ("~ lambda== 2.5 ~")"),
cex = 0.9)
# degraus e saltos
segments(xis[-10], acumulada[-10], xis[-1], acumulada[-10], lwd = 2)
segments(xis[-10], c(0, acumulada[-c(9,10)]),
xis[-10], acumulada[-10],
lty = 2)
segments(-1, 0, 0, 0, lwd = 2)
# pontos
points(xis, c(0, acumulada[-10]), pch = 21, cex = 0.7, bg = "white")
points(xis, acumulada, pch = 20, xpd = TRUE)
Observe que \(f_X(x_i) = \operatorname{P}\{X = x_i \}\) é igual ao salto que a função de distribuição acumulada \(F_X\) dá no ponto \(x_i\); por exemplo, \[ f_X(2)= \operatorname{P}\{X = 2 \} = 0.257 = F_X(2) - F_X(2-) \] ou seja,
ppois(2, taxa) - ppois(1, taxa)## [1] 0.2565156valor que também é obtido pelo comando dpois(2, taxa). Verifique.
De modo geral, \(f_X(x_i) = \operatorname{P}\{X = x_i \} = F_X(x_i) - F_X(x_i-)\), onde \(F_X(a-) = \lim_{x \to a-}F_X(x)\).
A função qpois(p, lambda)fornece o valor da função de distribuição acumulada \(F_X^{-1}(p), ~ 0\leq p \leq 1\) de uma distribuição de Poisson com taxa \(\lambda\). Seja uma variável aleatória \(X\), com distribuição de Poisson com parâmetro \(\lambda = 2,5\). Sua mediana é \(\nu = F_X^{-1}(0.5) = 2\), pois:
qpois(0.5, 2.5)## [1] 2que está contida no intervalo \(\lambda -\ln 2\leq \nu <\lambda + {\frac {1}{3}}\), ou seja, entre (1.807, 2.833).
Exemplos:
Ex1 - Suponha que uma bola de tênis, da quadra ao lado, caia em seu quintal a uma taxa média de 3 bolas por hora durante o dia. Qual é a probabilidade de que 10 ou menos bolas de tênis caiam em seu quintal durante a tarde, supondo que a tarde tenha 5 horas de duração?
# media de bolas perdidas por tarde
(media <- 3 * 5)## [1] 15 # probabilidade de 10 bolas em uma tarde típica
ppois(10, media)## [1] 0.1184644A média é \(3 \text{ bolas/h} \times 5 \text{ h} = 15\) bolas durante toda a tarde. A probabilidade solicitada é 0.1185.
Ex2 - Se, em sua vizinhança, um alarme dispara durante a madrugada taxa média de 1 a cada 4 horas, qual é a probabilidade de que pelo menos um dispare na próxima hora (questão escrita durante uma madrugada)?
# média de disparos (de alarme) por hora
(media <- 1/4)## [1] 0.25 # probabilidade de pelo menos um dispare (1 - P(nenhum dispara))
1 - ppois(0, 0.25)## [1] 0.2211992