Distribuições de Probabilidade no R

O R inclui algumas operações com as distribuições de probabilidade. No caso, existe 4 operações básica indicada pela seguinte letras:

Distribuição Binomial

  • Seja \(X\) o número de sucessos obtidos na realização de \(n\) ensaios de Bernoulli independentes

\[ \mathbb{P}[X=k]=\left(\begin{array}{c}n\\k\end{array}\right) p^k(1-p)^{n-k}. ~~~~~~~~~ k = 0,1,\ldots,n \]

Uma moeda honesta é lançada quatro vezes. Assim, o nosso espaço amostral \((\Omega)\) pode ser apresentado pela tabela a seguir:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
coroa coroa coroa coroa cara coroa coroa cara cara cara cara coroa cara cara cara cara
coroa coroa coroa cara coroa coroa cara cara coroa coroa cara cara cara coroa cara cara
coroa coroa cara coroa coroa cara cara coroa coroa cara coroa cara cara cara coroa cara
coroa cara coroa coroa coroa cara coroa coroa cara coroa coroa cara coroa cara cara cara

  • Poderíamos está interessado em responder a seguinte pergunta: Ao lançarmos uma moeda quatro vezes, qual a probabilidade de ocorrer cara, exatamente, duas vezes.

  • Para responder a este questionamento, precisamos montar a seguinte trabela de frequências:

Quantidade de Caras Frequência absoluta Frequência Relativa Frequência Relativa Acumulada
zero 1 0,0625 0,0625
uma 4 0,2500 0,3125
duas 6 0,3750 0,6875
três 4 0,2500 0,9375
quatro 1 0,0625 1,0000
Total 16 1,0000

Nota-se por meio da tabela de frequências que podemos dizer que a probabilidade de ocorrer exatamente duas caras em 4 lançamentos de uma moeda é de 0,375.

Outra conclussão: A probabilidade de ocorrer no máximo três caras em quatro lançamentos de uma moeda é de 0,9375.

moeda <- data.frame('cara'=
             rep(c('0',"1","2","3","4"), 
                  c(1, 4, 6, 4, 1)))
require(ggplot2)
## Loading required package: ggplot2
ggplot(moeda, aes(x = cara, y=..count../sum(..count..), fill = cara)) +
  geom_bar(width=0.1) +
  labs(title = "Lançamento de uma moeda quatro vezes",
     x = "Quantidade de caras",
     y = "Probabilidade",
     fill='QTD caras')

  • Note que dar bastante trabalho montar uma tabela de frequência. A medida que aumenta a quantidade de lançamentos de uma moeda, o espaço amostral cresce significativamente.

  • Para resolver este e outros eventuais problemas, foi construída a distribuição de probabilidade Binomial. A qual deve ser utilizada que estamos diante de um experimente em que temos apenas dois resultados possíveis: sucesso e fracasso.

No r, vamos apresentar como calcular a probabilidade de:

  1. Ocorrer exatamente 2 caras em 4 lançamentos de uma moeda.
dbinom(x=3, #Calcula a probabilidade de P(X=x)
       size=4, #Quantidade total de lançamentos
       prob=0.5, #Probabilidade inicial de ocorrer o sucesso
       log = FALSE)
## [1] 0.25
  1. Ocorrer no máximo 3 caras em 4 lançamentos de uma moeda.
 pbinom(q=3, #Quantidade de caras
        size=4, #Quantidade total de lançamentos
        prob=0.5, #Probabilidade inicial de ocorrência (lançar uma moeda uma única vez, qual a probabilidade de ocorrer cara)
        lower.tail = TRUE #P[X<= x]
        )
## [1] 0.9375
  1. Caso seja necessário, podemos utilizar a função plot para construir o gráficos de probabilidades
caras <- 0:4
probabilidade <- dbinom(x=caras,   # Quantidade de sucessos
                        size = 4,  # Quantidade de lançamentos
                        prob=0.5)  # Probabilidade a priori de sucesso

probabilidade
## [1] 0.0625 0.2500 0.3750 0.2500 0.0625
plot(caras, probabilidade, 
     type='h',   # Desenha uma linha vertical
     col='red',   # Cor da linha
     lwd=3)      # Espessura da linha/ponto

  • Exemplo : Suponha que numa linha de produção a probabilidade de se obter uma peça defeituosa (sucesso) é \(p=0,5\). Toma-se uma amostra de 10 peças para serem inspecionadas. Em que X é o número de peças defeituosa produzidas em único dia.

  • Qual é a probabilidade de um remessa com 10 peças conter pelo menos quatro peças defeituosa.

  • Note que, Pelos quatro significa, \(P[X \geq 4]\), isto é equivalente a \(1 - P[X<4]\) ou \(1 - P[X \leq 3]\).

Vamos ao R:

\[P[X \geq 4]\]

 pbinom(q=3, #Quantidade de peças defeituosas
        size=10, #Quantidade total de peças
        prob=0.5, #Probabilidade inicial de peça defeituosa
        lower.tail = FALSE #P[X> x]
        )
## [1] 0.828125

Resolva

Considere nascimentos de 4 filhotes de coelhos de uma determinada raça. Nesta raça há um distúrbio genético e a probabilidade de nascer fêmea é 5/8. Sendo X a ocorrência de fêmeas e utilizando a distribuição binomial obter:

  1. Construa um gráfico de probabilidades
femeas <- 0:4
probabilidade <- dbinom(x=femeas,   # Quantidade de sucessos
                        size = 4,   # Quantidade de nascimento
                        prob=5/8)   # Probabilidade a priori de sucesso

probabilidade
## [1] 0.01977539 0.13183594 0.32958984 0.36621094 0.15258789
plot(femeas, probabilidade, 
     type='h',   # Desenha uma linha vertical
     col='red',   # Cor da linha
     lwd=3)      # Espessura da linha/ponto

  1. Qual a probabilidade de nascer pelo menos três fêmeas. Dica: \(P[X \geq 3] = P[X > 2]\)
 pbinom(q=2,         # Quantidade de fêmeas
        size=4,      # Quantidade total de filhores
        prob=5/8,    # Probabilidade inicial de fêmea
        lower.tail = FALSE #P[X> x]
        )
## [1] 0.5187988

Distribuição Poisson

  • Seja \(X\) uma variável aleatória discreta que representa o número de sucesso. \[\mathbb{P}(X=k)=\frac{e^{-\lambda}\lambda^k}{k!}.\]

Problema: Nos últimos 30 dias foram registrados a quantidade de acidentes por dia em uma determinada cidade.

acidente <- data.frame('acd'= c(9, 6, 9, 11, 10, 10, 6, 10, 9, 4, 8, 10,
                                10, 7, 9, 11, 4, 6, 11, 8, 5, 3, 5, 9, 6))

table(acidente$acd)
## 
##  3  4  5  6  7  8  9 10 11 
##  1  2  2  4  1  2  5  5  3
ggplot(acidente, aes(x = as.factor(acd), y=..count../sum(..count..),
                     fill = as.factor(acd))) +
  geom_bar(width=0.1) +
  labs(title = "Acidentes por dia",
     x = "Quantidade de acidentes",
     y = "Frequência",
     fill='QTD acidentes')

A questão aqui é que não temos condições de afirmar qual a quantidade máxima de acidentes que podem ocorrer em um único dia.

  • Desta forma, o modelos de Poisson é o sugerido para dados de contagem dentro de um intervalo (temporal ou espacial) fixo.

Olhando para os dados, temos que a média e o desvio padrão da quantidade de acidentes por dia são de:

mean(acidente$acd)
## [1] 7.84
sd(acidente$acd)
## [1] 2.44404

  • Observando o gráfico parece razoável aformar que: a probabilidade de ocorrer 9 acidentes em um único dia é de 0,20.

  • Mas, olhando para a amostra, não somos capazes de afirmar qual a probabilidade de ocorrer 12 acidentes em um único dia.

  • Para responder a este questionamente, temos que fazer o uso de um modelo probabilístico, neste caso, Poisson.

  • Neste modelo, uma estimativa relacionada ao parâmetro \(\lambda\) deve ser a média amostral. Logo, podemos utilizar \(\lambda=8\).

Questionamento: qual a probabilidade de ocorrer 12 acidentes em um único dia.

dpois(x = 12,
      lambda = 8)
## [1] 0.0481268
acid_dia <- 0:20
probs <- dpois(x = acid_dia,
              lambda = 8)

probs
##  [1] 0.0003354626 0.0026837010 0.0107348041 0.0286261442 0.0572522885
##  [6] 0.0916036616 0.1221382155 0.1395865320 0.1395865320 0.1240769173
## [11] 0.0992615338 0.0721902064 0.0481268043 0.0296164949 0.0169237114
## [16] 0.0090259794 0.0045129897 0.0021237599 0.0009438933 0.0003974287
## [21] 0.0001589715
plot(acid_dia, probs, 
     type='h',   # Desenha uma linha vertical
     col='red',   # Cor da linha
     lwd=4)      # Espessura da linha/ponto

Exemplo: Uma empresa de telefonia recebe, em média, cinco chamadas por minuto. Supondo que a distribuição de Poisson seja adequada nessa situação, obter:

  1. Gráfico da distribuição de probabilidade
tel_min <- 0:10
probs <- dpois(x = tel_min,
              lambda = 5)

probs
##  [1] 0.006737947 0.033689735 0.084224337 0.140373896 0.175467370
##  [6] 0.175467370 0.146222808 0.104444863 0.065278039 0.036265577
## [11] 0.018132789
plot(tel_min, probs, 
     type='h',   # Desenha uma linha vertical
     col='red',   # Cor da linha
     lwd=4)      # Espessura da linha/ponto

  1. A probabilidade de que a empresa não receba chamada durante um intervalo de um minuto.
dpois(x = 0, lambda = 5)
## [1] 0.006737947

No caso, para outras distribuição temos que:

# Binomial Negativa        dnbinom(n,size, prob)

# Geométrica               dgeom(n,prob)

# Hipergeometrica          dhyper(nn, m, n, k)

Distribuição Normal

Uma variável aleatória contínua \(X\) tem distribuição Normal se sua função densidade de probabilidade for dada por:

\[f(x)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left[-\frac{1}{2}\left(\frac{x-\mu}{\sigma}\right)^2\right], \quad x\in(-\infty,\infty).\]

A seguir temos dados de variáveis climatológicas.

Link para download: https://www.dropbox.com/s/zx6311dpfcso2c0/clim_dou.csv?dl=0

clim_dou <- read.csv2("C:/Users/Carol/Dropbox/UFGD/2019.01_Disciplinas/Topicos de Estatistica/4_Aula/clim_dou.csv")

head(clim_dou)
##         data Tem_max Tem_min Chuva  Ano Meses
## 1 01/01/1980    25.0    11.0     0 1980     1
## 2 02/01/1980    28.2    12.1     0 1980     1
## 3 03/01/1980    30.5    14.8     0 1980     1
## 4 04/01/1980    31.9    18.7     0 1980     1
## 5 05/01/1980    31.8    18.1     0 1980     1
## 6 06/01/1980    30.2    18.2     0 1980     1

Vamos selecionar apenas o mês de janeiro

jan = subset(clim_dou, Meses==1)
  • Vamos construir o gráfico de histograma para a temperatura máxima diária
ggplot(jan, aes(x = Tem_max, ..density..)) +
  
  geom_histogram(binwidth = 1,     # Amplitude da classe
                 fill = 'dodgerblue',
                 color = 'black') +
  
  labs(title = "Distribuição da temperatura máxima diária",
       x = "Temperatura Máxima (º C)",
       y = "Densidade") +
  
  geom_density(alpha = 0.5)+     # Linha de densidade
  
  theme_light()

Problema: Como calcular a probabilidade, em um único dia, da temperatura máxima ser superior a 36ºC.

  • Neste problema temos o cálculo de área, o qual envolve o cálculo de integral. Mas qual função devo integrar? Para responder a este questionamento, vamos verificar se a função de densidade Normal é adequada para modelar esta variável.
ggplot(jan, aes(x = Tem_max)) +
  
  geom_histogram(aes(y = ..density..),
                 binwidth = 1,     # Amplitude da classe
                 fill = 'dodgerblue',
                 color = 'black') +
  
  labs(title = "Distribuição da temperatura máxima diária",
       x = "Temperatura Máxima (º C)",
       y = "Densidade") +
  
  geom_density(alpha = 0.5)+     # Linha de densidade
  
  stat_function(fun = dnorm, color='red', size = 2,
                args = list(mean = mean(jan$Tem_max),
                            sd = sd(jan$Tem_max)))+   # Linha de densidade da distribuição Normal
  
  theme_light()

Então, o cálculo de probabilidade se resumo a resolver a integral:

mean(jan$Tem_max)
## [1] 31.7134
sd(jan$Tem_max)
## [1] 2.412448
  • Qual a probabilidade, em um único dia, da temperatura máxima ser superior a 36ºC:

\[P[X \geq 36] = \int_{36}^{\infty}\frac{1}{\sqrt{2\pi 2,41}}\exp\left[-\frac{1}{2}\left(\frac{x-31,71}{2,41}\right)^2\right]dx \]
Vamos ao R

pnorm(q = 36,        # x=36
      mean = 31.71,  # média
      sd = 2.41,     # Desvio padrão
      lower.tail = FALSE # Calcula P[X > x]
      )
## [1] 0.03753119

Assim, a temperatura máxima em um único dia exceder 36 ºC é de 0,04

  • Calcule a probabilidade da temperatura máximia em um único dia em janeiro ser:
  1. Inferior a 30 ºC
pnorm(q = 30,        # x=30
      mean = 31.71,  # média
      sd = 2.41,     # Desvio padrão
      lower.tail = TRUE # Calcula P[X ??? x]
      )
## [1] 0.2389936
  1. Está entre 30 ºC e 35 ºC.
pnorm(q = c(30, 35),        # x=30
      mean = 31.71,  # média
      sd = 2.41,     # Desvio padrão
      lower.tail = TRUE # Calcula P[X ??? x]
      )
## [1] 0.2389936 0.9138963
0.9138963 - 0.2389936
## [1] 0.6749027

Resolva: Para a temperatura mínima diária de janeiro faça:

1.) Construa um gráfico de histograma em que seja possível analisar se a distribuição Normal é adequada para modelar a temperatura mínima.

ggplot(jan, aes(x = Tem_min)) +
  
  geom_histogram(aes(y = ..density..),
                 binwidth = 1,     # Amplitude da classe
                 fill = 'green',
                 color = 'gray60') +
  
  labs(title = "Distribuição da temperatura mínima diária",
       x = "Temperatura mínima (º C)",
       y = "Densidade") +
  
  geom_density(alpha = 0.5)+     # Linha de densidade
  
  stat_function(fun = dnorm, color='red', size = 2,
                args = list(mean = mean(jan$Tem_min),
                            sd = sd(jan$Tem_min)))+   # Linha de densidade da distribuição Normal
  
  theme_light()

  1. Calcule a probabilidade da temperatura mínima em um único dia de janeiro ser:
  1. Inferior a 20 ºC
mean(jan$Tem_min)
## [1] 21.01679
sd(jan$Tem_min)
## [1] 1.60775
pnorm(q = 20,        # x=20
      mean = 21.02,  # média
      sd = 1.61,     # Desvio padrão
      lower.tail = TRUE # Calcula P[X ??? x]
      )
## [1] 0.2631904
  1. Superior a 22 ºC
pnorm(q = 22,            # x=22
      mean = 21.02,      # média
      sd = 1.61,         # Desvio padrão
      lower.tail = FALSE # Calcula P[X > x]
      )
## [1] 0.2713631

Distribuição Exponencial

A variável aleatória \(X\) tem distribuição Exponencial com parâmetro \(\lambda\),se tiver função densidade de probabilidade dada por:

\[f(x)=\left\{\begin{array}{l}\lambda e^{-\lambda x} \ \hbox{se} \ x\geq 0\\0 \ \hbox{se} \ x \ < \ 0\end{array}\right.\]

Exemplo: O tempo até a falha do ventilador de motores a diesel tem uma distribuição Exponencial com parâmetro \(\lambda = \frac{1}{28700}\) horas. Qual a probabilidade de um destes ventiladores falhar nas primeiras 24000 horas de funcionamento?

x<- 24000
razao<- 1/28700
prob<-pexp(x,razao);prob
## [1] 0.5666619

Ou seja, a probabilidade de um destes ventiladores falhar nas primeiras $ 24000$ horas de funcionamento é de, aproximadamente, 56,7%

Quantil da distribuição

qf(0.975,12,10)
## [1] 3.620945
qt(0.9,10)
## [1] 1.372184
qchisq(0.975,12)
## [1] 23.33666
# Exemplo da normal : Qual deve ser o tempo de prova de modo a permitir que 95% dos vestibulandos terminem no prazo estipulado?
qnorm(0.95,120,15)
## [1] 144.6728

Gerar número aleatório de uma distribuição de interesse

n=5
t <- rt(n,10);t
## [1] -0.01347506 -1.15017030 -0.44663557  1.25372468  0.38960113
normalp <- rnorm(10);normalp
##  [1] -0.52916230 -2.48241327  1.08597494 -0.33566988  0.09418351
##  [6]  0.57853408 -0.56159653  1.64843190 -1.00619723  1.19489573
normal <- rnorm(10,2,3); normal
##  [1]  2.3790584  0.4720149 -0.2871017 -3.4771221  5.5861796 -0.2056282
##  [7]  0.3092887 -1.5552943  4.9377231  4.4792388
bino<-rbinom(10,10,0.3);bino
##  [1] 2 4 2 3 3 2 3 2 4 4
#Uniforme          runif(n, min=0, max=1) 
#Exponencial       rexp(n, rate=1) 
#T- Student        rt(n, df) 'Student' (t)
#F                 rf(n, df1, df2) 
#Quiquadrado       rchisq(n, df) 
#Gama              rgamma(n, shape, scale=1) 
#Beta              rbeta(n, shape1, shape2) beta
#Log- Normal       rlnorm(n, meanlog=0, sdlog=1) lognormal
#Cauchy            rcauchy(n, location=0, scale=1) Cauchy