aula 3

aula de hoje: probabilidades

1. a função sample()

2. jogando dados

3. tirando bolas

4. decisão por pênaltis

5. aplicações simples do teorema de Bayes: vida extragaláctica

6. o problema de Monty Hall

7. o dilema do prisioneiro

8. combinação de probabilidades

9. exercícios

---

---

1 a função sample()

A função sample() permite amostrar aleatoriamente um vetor

Ela pode ser usada com substituição (replacement = TRUE) ou não (replacement = FALSE)

Imagine que você gera uma amostra de bolas coloridas e as coloca numa caixa: com substituição, você tira uma bola e a retorna para a caixa (pode então, haver repetição desse elemento);
sem substituição, não

# reprodutibilidade: 
# devemos estabelecer a semente do gerador de números aleatórios para termos
# resultados reprodutíveis
set.seed(321)

# vetor com números inteiros de 1 a 10
v = 1:10
v

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

# amostragem do vetor
# o default é replacement = FALSE
sample(v)

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

# de novo
sample(v)

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

# amostra com 5 elementos de uma amostra de inteiros de 1 a 10
sample(v, size = 5)

[1]  7  2  4  8 10

# outra amostra com 5 elementos de uma amostra de inteiros de 1 a 10
sample(v, size = 5)

[1]  3  7  1  8 10

# agora permitindo substituição (um elemento pode aparecer mais de uma vez)
sample(v, size = 15, replace = TRUE)

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

# de novo
sample(v, size = 15, replace = TRUE)

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

# de novo
sample(v, size = 6, replace = TRUE)

[1] 3 9 6 7 2 7

# tente rodar o comando:
# sample(v, size = 15, replace = FALSE)
# e entenda porque dá erro

---

---

2 jogando dados

Jogamos um dado de 6 faces: qual é a probabilidade do resultado ser par?

Vamos escrever o espaço amostral do problema, o espaço com todos os resultados possíveis: \[ E = \{ 1,2,3,4,5,6\} \] a dimensão do espaço amostral é 6

resultados pares: \(2,4,6\) (dimensão 3)

probabilidade do resultado ser par: \(3/6 = 1/2\)

vamos fazer uma simulação para verificar isso

vamos jogar o dado \(n=10000\) vezes e ver a frequência de resultados pares:

# reprodutibilidade: 
# devemos estabelecer a semente do gerador de números aleatórios para termos 
# resultados reprodutíveis
set.seed(123)

# vamos criar uma função para gerar uma amostra de n pontos de uma distribuição uniforme discreta entre 1 e k 
# amostragem com substituição
rdu<-function(n,k){sample(1:k,n,replace=T)}

# exemplos
rdu(10,5)

 [1] 3 3 2 2 3 5 4 1 2 3

rdu(5,10)

[1] 5 3 9 9 9

# vamos agora jogar dados
# número de simulações
n=10000

# simulação: jogando um dado n vezes
# para um dado de 6 faces, k=6
dados <-rdu(n,6)

# primeiros valores do vetor n-dimensional
head(dados)

[1] 5 3 2 2 1 6

# função para ver se um número é par (TRUE ou FALSE)
ehpar <- function(x){x %% 2 == 0}
# x %% y: x módulo y - dá o resto da divisão de um número por outro
# note o operador lógico ==

par = ehpar(dados)

# note que o vetor par contém o valor lógico (TRUE ou FALSE) dos números simulados: TRUE se for par e FALSE se for ímpar

# primeiros valores do vetor n-dimensional
head(par)

[1] FALSE FALSE  TRUE  TRUE FALSE  TRUE

# número de pares:
npar = length(par[par == 'TRUE'])

# fração de pares:
npar/n

[1] 0.4954

# a diferença em relação ao valor esperado (0.5) é devido às flutuações estatísticas da amostragem

---

---

3. tirando bolas

Uma caixa contém 4 bolas azuis, 2 vermelhas e 3 amarelas.

• Tiramos uma bola ao acaso e a retornamos à caixa.

• Repetimos isso 3 vezes.

• Qual é a probabilidade de termos tirado 2 bolas azuis e 1 amarela?

Vamos representar as bolas azuis com o número 1, as vermelhas com o 2 e as amarelas com o 3. Assim, a caixa pode ser representada pelo conjunto:

caixa <- c(1,1,1,1,2,2,3,3,3)

Temos 3 possibilidades de tirar 2 bolas azuis e uma amarela: (1,1,3), (1,3,1) e (3,1,1), cada uma com probabilidade \(4/9 \times 4/9 \times 3/9 = 48/729\). Logo, a probabilidade de tirarmos 2 bolas azuis e 1 amarela é \(3 \times 48/729 \simeq 0.1975\).

Vamos verificar este resultado fazendo \(n=10000\) simulações onde, de cada vez, tiramos 3 bolas:

# reprodutibilidade: 
# devemos estabelecer a semente do gerador de números aleatórios para termos resultados reprodutíveis
set.seed(123)

# número de simulações
n = 10000

# inicialização do contador
ncerto = 0
nc=0
# o comando for() cria um loop sobre o que está dentro dos {}
for (i in 1:n){
# amostragem com substituição: retorna-se a bola à caixa
# tira-se 3 bolas:
a = sample(caixa, size = 3, replace = TRUE) 
# verifica-se se o número de bolas azuis (1) é 2 e o de amarelas (3) é 1:
# usamos operadores lógicos == (igual) e & (AND)
if(length(a[a == 1]) == 2 & length(a[a == 3]) == 1) {ncerto = ncerto + 1}
}
# fração simulada:
ncerto/n

[1] 0.1953

repetindo:

# precisamos reinicializar a variável ncerto
ncerto = 0
for (i in 1:n){
a = sample(caixa, size = 3, replace = TRUE) 
if(length(a[a == 1]) == 2 & length(a[a == 3]) == 1) {ncerto = ncerto + 1}
}
# fração simulada:
ncerto/n

[1] 0.1981

Compare o resultado dessas duas simulações com o valor teórico. As diferenças mostram o efeito da amostragem.

Vamos agora considerar o mesmo problema onde a bola tirada não retorna à caixa. Nesse caso, a probabilidade de cada uma das 3 possibilidades fica

• \(P(1,1,3) = 4/9 \times 3/8 \times 3/7 = 36 / 504\)

• \(P(1,3,1) = 4/9 \times 3/8 \times 3/7 = 36 / 504\)

• \(P(3,1,1) = 3/9 \times 4/8 \times 3/7 = 36 / 504\)

e a probabilidade total é \(3 \times 36/504 \simeq 0.2142...\)

Vamos verificar isso com simulações e para ver como o resultado depende de \(n\), vamos usar function() para definir uma simulação com \(n\) amostras de 3 bolas:

simula <- function(n){
ncerto = 0
for (i in 1:n){
# sem replacement
a = sample(caixa, size = 3, replace = FALSE) 
if(length(a[a == 1]) == 2 & length(a[a == 3]) == 1) {ncerto = ncerto + 1}
}
# fração simulada:
return(ncerto/n)
}

a função simula(n) retorna a fração com que se obtém 2 bolas azuis e uma amarela em n simulações.

exemplos:

simula(10)

[1] 0.2

simula(100)

[1] 0.16

simula(1000)

[1] 0.202

simula(1000000)

[1] 0.214503

Note como a fração estimada se aproxima da esperada conforme o número de simulações aumenta.

---

---

4. decisão por pênaltis

Em um jogo de futebol A e B vão chutar um pênalti cada um. A taxa de acerto deles é 0.7 e 0.8, respectivamente.

qual é a probabilidade de que pelo menos um deles faça gol?

temos que \(P(A)=0.7\) e \(P(B)=0.8\)

a probabilidade de que pelo menos um deles faça gol é 1 menos a probabilidade de que nenhum faça gol: \[ 1-[1-P(A)] \times [1-P(B)]=1-0.06=0.94\]

• qual é a probabilidade de que só B faça gol? \[ (1-P(A))P(B)=0.3\times 0.8=0.24 \]

• qual é a probabilidade de que só um dos jogadores faça gol? \[ P(A)(1-P(B))+P(B)(1-P(A))=0.7\times 0.2+0.8\times 0.3=0.38\]

Vamos testar esses resultados com uma simulação. Como podemos fazer isso?

A ideia básica do método de Monte Carlo aplicado a este problema é a seguinte: sorteie um número uniformemente distribuído entre 0 e 1, r. Se \(P(A) > r\), consideramos que A fez um gol, em caso contrário, não. O mesmo procedimento vale para B.

# reprodutibilidade
set.seed(123)

# número de simulações: em cada uma, A e B dão um chute cada um
nsim=1000000

# probabilidade de marcar o gol
pa=0.7
pb=0.8

# inicialização
npm1=0 # numero de simulações onde pelo menos um fez gol
nsob=0 # numero de simulações onde só B fez gol
nso1=0 # numero de simulações onde só um dos jogadores fez gol

# simulacões
# vamos definir duas variáveis, A e B, para representar se A ou B marcam
for(i in 1:nsim){

# A chuta
A=0
r = runif(1)
if(pa > r)A=1

# B chuta
B=0
r = runif(1)
if(pb > r)B=1

# note os operadores lógicos OR: | e AND: &
if(A == 1 | B == 1)npm1=npm1+1
if(A + B == 1)nso1=nso1+1
if(A == 0 & B == 1)nsob=nsob+1
}

# resultado
c(npm1/nsim,nso1/nsim,nsob/nsim)

[1] 0.939858 0.379649 0.239819

# ou, arredondando,
c(round(npm1/nsim,4),round(nso1/nsim,3),round(nsob/nsim,3))

[1] 0.9399 0.3800 0.2400

compare com os valores obtidos analiticamente: 0.94, 0.38 e 0.24

---

---

5. aplicações simples do teorema de Bayes: vida extragaláctica

Em uma galáxia muito, muito distante:
\(-\) 5% dos planetas têm água líquida: P(agua) = 0.05
\(-\) 1% dos planetas têm alguma forma de vida: P(vida) = 0.01
\(-\) A probabilidade de haver vida dado que o planeta tem água líquida é 10%: P(vida∣agua) = 0.10

Qual é a probabilidade de um planeta ter água líquida, dado que ele tem vida? P(agua∣vida) = ?

Solução analítica (Bayes): \[P(agua∣vida) = {P(vida∣agua) P(agua) \over P(vida)}\]

Substituindo os valores: \[P(agua∣vida) = {0.10 \times 0.05 \over 0.01} = 0.5\] Nessa galáxia, metade dos planetas com vida tem água líquida.

simulação: vamos
1. Criar uma população grande de planetas.
2. Sortear quais têm água.
3. Sortear quais têm vida condicionado a terem ou não água.
4. Medir empiricamente P(água∣vida)

Para a simulação ser consistente, precisamos conhecer a probabilidade de vida se não tiver água no planeta, P(vida|!água).
Para isso usamos a lei da probabilidade total: \[ P(vida)=P(vida∣agua)P(agua)+P(vida∣!agua)P(!agua)\] e, portanto, \[ P(vida∣!agua) = {P(vida)−P(vida∣agua)P(agua)​ \over 1-P(agua)} = {0.01-0.10 \times 0.05 \over 0.95} = {0.005 \over 0.95} \simeq 0.00526\] Ou seja: nessa galáxia, vida sem água é raríssima.

# reprodutibilidade
set.seed(42)

N <- 1e6         # Número de planetas simulados

# -------------------------------
# 1) Probabilidades dadas
# -------------------------------
p_agua <- 0.05
p_vida <- 0.01
p_vida_com_agua <- 0.10

# -------------------------------
# 2) Probabilidade derivada:
#    P(vida | não água)
# -------------------------------
p_vida_sem_agua <-
  (p_vida - p_vida_com_agua * p_agua) / (1 - p_agua)

p_vida_sem_agua

[1] 0.005263158

# -------------------------------
# 3) Sorteio de água líquida
# -------------------------------
agua <- rbinom(N, size = 1, prob = p_agua) == 1

# -------------------------------
# 4) Sorteio de vida condicionado à água
# -------------------------------
vida <- logical(N)

vida[agua] <-
  rbinom(sum(agua), size = 1, prob = p_vida_com_agua) == 1

vida[!agua] <-
  rbinom(sum(!agua), size = 1, prob = p_vida_sem_agua) == 1

# -------------------------------
# 5) Verificações
# -------------------------------
cat("P(água) simulado =", mean(agua), "\n")

P(água) simulado = 0.04997 

cat("P(vida) simulado =", mean(vida), "\n")

P(vida) simulado = 0.010263 

cat("P(vida|água) simulado =", mean(vida[agua]), "\n")

P(vida|água) simulado = 0.1042225 

# -------------------------------
# 6) Estimar P(água | vida)
# -------------------------------
p_agua_com_vida_sim <- mean(agua[vida])
cat("P(água | vida) simulado =", p_agua_com_vida_sim, "\n")

P(água | vida) simulado = 0.507454 

# -------------------------------
# 7) Valor teórico (Bayes)
# -------------------------------
p_agua_com_vida_teo <-
  (p_vida_com_agua * p_agua) / p_vida

cat("P(água | vida) teórico =", p_agua_com_vida_teo, "\n")

P(água | vida) teórico = 0.5 

---

---

6. o problema de Monty Hall

você tem 3 portas na sua frente: uma contém uma Ferrari e as duas outras um bode em cada uma

MH pede para você escolher uma delas: você escolhe, por exemplo, a 1

antes de abrir a porta que você escolheu, MH (que sabe onde o carro está), escolhe uma das portas com um bode e a abre; sobra assim uma outra porta fechada

MH então te pergunta: quer trocar de porta ou não quer? Você quer trocar a porta que escolheu por esta que sobrou ou não?

o que é mais vantajoso fazer?

vamos fazer uma simulação do problema:

(baseado em https://www.r-bloggers.com/2020/10/monty-hall-a-programmers-explanation/)

# reprodutibilidade: 
# devemos estabelecer a semente do gerador de números aleatórios para termos resultados reprodutíveis
set.seed(123)

# número de "jogadas"
n = 10000

# vamos chamar as portas de '1', '2','3' 
portas = c('1', '2', '3')

# o vetor carro contém a porta onde o carro está em cada simulação
carro =  sample(portas, size = n, replace = TRUE)

# vamos ver os primeiros valores:
head(carro)

[1] "3" "3" "3" "2" "3" "2"

# vamos simular a porta que MH abre e vamos supor, sem perda de generalidade, que você abriu a '1'
# se o carro está em '1', Monty vai abrir ou a '2' ou a '3', ao acaso; se não, vai abrir a porta que tem o bode
# comando ifelse(test, yes, no)

porta <- ifelse(
  carro == '1',                                  # se o carro está em 1
  sample(c('2', '3'), size = n, replace = TRUE), # y: escolhe ao acaso 
  ifelse(carro == '2', '3', '2'))                # n: escolhe a porta com o bode

# vamos ver a probabilidade de se ganhar se você não trocar de porta, que é dada pela fração de vezes em que o carro está atrás da porta '1':
sum(carro == '1') / n

[1] 0.3335

# como esperado é mais ou menos 1/3

# vamos ver a probabilidade de se ganhar se você trocar de porta, que é a fração de vezes em que a porta fechada tem o carro:
porta_fechada <- ifelse(porta == '2', '3', '2')
sum(carro == porta_fechada) / n

[1] 0.6665

# o carro ou está atrás da porta fechada ou da porta '1'; você ganha trocando só se o carro não estiver em '1'
# assim, a última expressão é equivalente a
sum(carro != '1') / n

[1] 0.6665

# que dá o mesmo resultado

a probabilidade de ganhar o carro sem trocar de porta é a probabilidade de que o carro esteja atrás da porta ‘1’, 1/3, e a probabilidade de se ganhar trocando de porta é a probabilidade de que o carro não esteja em ‘1’, que é 2/3

---

---

7. o dilema do prisioneiro

Um carrasco chama um prisioneiro e lhe mostra 4 portas: três conduzem à forca e uma à liberdade. O carrasco pede ao prisioneiro para escolher uma porta mas, antes de abri-la, abre uma outra, que conduz à forca, e pergunta se o prisioneiro quer ficar com a porta que escolheu ou se quer trocar por uma das outras duas que permaneceram fechadas. O que o prisioneiro deve fazer? Qual é a probabilidade dele ganhar a liberdade sem e com trocar de porta?

(este problema é análogo ao de Monty Hall)

---

---

8 combinação de distribuições de probabilidades

Exemplo 1 Suponha que uma variável \(\phi\) seja uniformemente distribuida entre 0 e \(2\pi\).
Qual é a distribuição esperada de \(\sin \phi\)?

Temos que \(P(\phi)\) é uniforme: \(P(\phi) = c\), onde \(c\) é uma constante

\(P(\phi)\) deve ser normalizada: \[ \int_0^{2 \pi} P(x)dx = 2 \pi c = 1\] logo, \(c= 1/(2 \pi)\) e \[ P(\phi) = \frac{1}{2 \pi}~~~~ (0 < \phi < 2\pi) \]

Seja \(y = \sin \phi\). Temos \(P(\phi)\) e queremos \(P(y)\)

Como \[ P(y)dy = P(\phi)d\phi ~\rightarrow ~ P(y) = P(\phi) \left| {d\phi \over dy} \right| \] e \[ \frac{dy}{d\phi} =\cos \phi =\sqrt{1-\sin^2 \phi} = \sqrt{1-y^2}, \] vem que \[ P(y) = \frac{1}{2 \pi \sqrt{1-y^2}} \] (note que \(-1 \le y \le 1\))

Vamos fazer uma simulação para verificar o resultado.

# reprodutibilidade
set.seed(321123)

# tamanho da amostra a ser simulada
n=100000L  #esse L significa que n é inteiro

# simulação de n valores uniformemente distribuídos entre min e max
fi=runif(n, min = 0, max = 2*pi)

# variável y
y=sin(fi)

#distribuição de y
h=hist(y)

# conteúdo da variável h:
h

$breaks
 [1] -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1  0.0  0.1  0.2  0.3  0.4
[16]  0.5  0.6  0.7  0.8  0.9  1.0

$counts
 [1] 14310  6220  4853  4252  3862  3558  3439  3354  3193  3184  3156  3166
[13]  3252  3415  3566  3813  3994  4864  6075 14474

$density
 [1] 1.4310 0.6220 0.4853 0.4252 0.3862 0.3558 0.3439 0.3354 0.3193 0.3184
[11] 0.3156 0.3166 0.3252 0.3415 0.3566 0.3813 0.3994 0.4864 0.6075 1.4474

$mids
 [1] -0.95 -0.85 -0.75 -0.65 -0.55 -0.45 -0.35 -0.25 -0.15 -0.05  0.05  0.15
[13]  0.25  0.35  0.45  0.55  0.65  0.75  0.85  0.95

$xname
[1] "y"

$equidist
[1] TRUE

attr(,"class")
[1] "histogram"

# para normalizar P(y)
hmin=min(h$counts) 

# gráfico
y=seq(-0.99,0.99,0.01)
yc=hmin/sqrt(1-y^2)
lines(y,yc,col='red',lwd=3)

---

---

---

Exemplo 2 Vamos agora considerar a distribuição da razão de duas variáveis gaussianas: se \(x\) e \(y\) são variáveis distribuidas como gaussianas \(N(0,1)\), qual é a distribuição de \(z=x/y\)?

Vamos ver a cara da distribuição de \(z\) fazendo uma simulação:

# reprodutibilidade
set.seed(456)

# tamanho da amostra
m=1000

# para simular números distribuídos como uma gaussiana usa-se a função rnorm()
# sintaxe: rnorm(m,média,desvio padrão)    o default é média 0 e desvio padrão 1
# simulação dos valores de x, y e z
x=rnorm(m)
y=rnorm(m)
z=x/y

# algumas estatísticas de z
summary(z)

      Min.    1st Qu.     Median       Mean    3rd Qu.       Max. 
-336.18158   -0.98541   -0.03031   -0.21466    0.94955  130.76489 

# gráfico: 2 painéis na mesma linha
par(mfrow=c(1,2))
hist(z,main="distribuição de z=x/y",ylab='frequência')

# vamos dar um zoom para -5 < z < 5
zz = z[z > -5 & z < 5]
hist(zz,main="zoom da distribuição de z",xlim=c(-5,5),breaks=100,xlab='z',ylab='frequência')

Pode-se mostrar que \(z\) obedece a uma distribuição de Cauchy: \(h(z) = {1 \over \pi} {1 \over 1+z^2}\)

# sequência de inteiros de 0 a 20 separados por 1
x = seq(-5,5,0.1) # numeros inteiros positivos
y = dcauchy(x) # distribuição de Cauchy com parâmetros 1 e 0
par(mfrow=c(2,2)) # panel de 2 x 2 figuras
plot(x,y,main="dist. de Cauchy",type="l",col='blue',ylim=c(0,0.4))

# distribuição gaussiana para comparação
lines(x,dnorm(x),col='red')
legend(-5, 0.4, legend=c("Gaussiana", "Cauchy"),
       col=c("red", "blue"), lty=1, cex=0.8,box.lty=0)
       
# distribuição cumulativa       
y = pcauchy(x) 
plot(x,pnorm(x),main="distr. cumulativa",type="l",col='red') # gaussiana
lines(x,y,col='blue') # cauchy

# quantis
x = seq(0,1,0.01)  # 100 números uniformemente distribuídos entre 0 e 1
y = qcauchy(x) 
plot(x,y,main="vetor de quantis",type="l",col='blue')
lines(x,qnorm(x),col='red')

# amostra aleatória
x=rcauchy(100)
hist(x,main="distr. aleatória ",col='blue',breaks=1000)

comparada à gaussiana, a distribuição de Cauchy tem “caudas pesadas”

---

---

9. exercícios da semana

• 1.Numa certa turma verifiquei que \(X\)% dos estudantes da classe disseram conhecer R e Python e \(Y\)% disseram conhecer Python. Qual é a probabilidade de um estudante que conhece Python conhecer também R? Inicialize a semente de números aleatórios do R com seu número USP e gere 1 número aleatório, \(g\): \(set.seed(SeuNumeroUsp); g = runif(1)\). Assuma que \(X=40+3 g\) e \(Y=60+5 g\).

• 2. Faça uma simulação para o slide:
     aplicações simples do teorema de Bayes: classificação visual

---

---

---

---