aula de hoje: probabilidades

combinando distribuições de probabilidades

simulações

exercícios: jogando dados

exercícios: tirando bolas

decisão por pênaltis

o problema de Monty Hall

o dilema do prisioneiro

exercícios

1 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.1816   -0.9854   -0.0303   -0.2147    0.9496  130.7649

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

# 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')

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”

2. Simulações

O R tem diversas funções que ajudam a simular distribuições de probabilidades.

Por exemplo, podemos gerar amostras de acordo com as distribuições de probabilidades mais comuns:

rnorm(): distribuição normal

runif(): distribuição uniforme

rpois(): distribuição de Poisson

Por exemplo, vamos gerar duas sequências independentes de 5 números com distribuição uniforme entre 0 e 10:

runif(5,min=0,max=10)

## [1] 8.3368227 2.9217463 4.9296249 0.5050883 9.6299453

runif(5,min=0,max=10)

## [1] 9.138157 5.513308 8.223068 6.410709 3.864730

# obtemos duas amostras diferentes de uma mesma população

mas usando a mesma semente obtém-se o mesmo resultado:

# mesma semente
set.seed(1234)
runif(5,min=0,max=10)

## [1] 1.137034 6.222994 6.092747 6.233794 8.609154

set.seed(1234)
runif(5,min=0,max=10)

## [1] 1.137034 6.222994 6.092747 6.233794 8.609154

# com outra semente, tem-se outra amostra
set.seed(1)
runif(5,min=0,max=10)

## [1] 2.655087 3.721239 5.728534 9.082078 2.016819

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

# 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]  7  2  3  8  1  5  6  9 10  4

# de novo
sample(v)

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

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

## [1] 5 2 9 6 1

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

## [1] 1 4 3 6 2

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

##  [1] 10  6  4  4 10  9  7  6  9  8  9  7  8  6 10

# de novo
sample(1:10, size = 15, replace = TRUE)

##  [1]  7  3 10  6  8  2  2  6  6  1  3  3  8  6  7

# de novo
sample(1:10, size = 6, replace = TRUE)

## [1] 6 8 7 1 4 8

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

3. exercícios: 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 n pontos com 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

4. exercícios: 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.

5. 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

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)

aula2

lsj

March 7, 2025

aula de hoje: probabilidades

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

2. Simulações

3. exercícios: jogando dados

4. exercícios: tirando bolas

5. decisão por pênaltis

6. o problema de Monty Hall

7. o dilema do prisioneiro

exercícios da semana