aula de hoje: probabilidades

distribuição gaussiana
simulações
visualização de distribuições bidimensionais
exercícios/simulações

1. DISTRIBUIÇÃO GAUSSIANA

distribuição gaussiana de média \(\mu\) e desvio padrão \(\sigma\): \[ P(x) = N(\mu,\sigma) = {1 \over \sqrt{2 \pi} \sigma} \exp\Big[-{(x-\mu)^2 \over 2 \sigma^2} \Big] \]

distribuição gaussiana em R:

dnorm(x, mean = 0, sd = 1, log = FALSE)

pnorm(q, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)

qnorm(p, mean = 0, sd = 1, lower.tail = TRUE, log.p = FALSE)

rnorm(n, mean = 0, sd = 1)

onde:

d- densidade: a fdp

p- probabilidade: a função de distribuição cumulativa

q- quantis: inverso da função de distribuição cumulativa

r- random: gera variáveis aleatórias com essa distribuição

densidade de probabilidades

x = seq(-4, 4, by = .1)
plot(x,dnorm(x,mean=0,sd=1),main="fdp gaussiana",type="l",col='blue')

distribuição cumulativa

a função pnorm(x) retorna a função de distribuição cumulativa da gaussiana, \(\int_{-\infty}^x P(x^\prime) d x^\prime\)

plot(x,pnorm(x,mean=0,sd=1),main="distribuição cumulativa da gaussiana",type="l")

percentil da distribuição cumulativa

qnorm é o inverso de pnorm e fornece o falor de x correspondente ao percentil p da distribuição normal

p = seq(0, 1, by = .01)
x = qnorm(p)
plot(p,x,main="qnorm",type="l")

geração de números aleatórios

rnorm gera um vetor com números distribuídos como uma gaussiana

# 100 números aleatórios gerados a partir de uma distribuição gaussiana de média 10 e desvio padrão 5
x=rnorm(100,mean=10,sd=5)
hist(x)

# media dos dados simulados
mean(x)

## [1] 10.42614

# desvio padrão dos dados simulados
sd(x)

## [1] 4.891585

um exemplo de modelagem gaussiana:

A massa do buraco negro no centro da Via Láctea é de 4.3 \(\pm\) 0.3 ( 1 sigma) milhões de massas solares. Qual é a probabilidade desta massa ser maior que 5 milhões de massas solares?

Supondo que a massa x (em unidades de milhões de massas solares) seja distribuída como P(x)=N(4.3,0.3), qual é a probabilidade de P(x>5)?

A distribuição cumulativa permite responder a esta questão:

x = seq(0, 10, by = .1)
plot(x,pnorm(x,mean=4.3,sd=0.3),main="distribuição cumulativa da pdf da massa do buraco negro",type="l")

probabilidade de P(x>5):

1-pnorm(5,mean=4.3,sd=0.3)

## [1] 0.009815329

A que massa do BN corresponde o percentil de 95%?

qnorm(0.95,mean=4.3,sd=0.3)

## [1] 4.793456

2. Simulações

O R tem diversas funções que ajudam em simulações.

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

rnorm(): distribuição normal
runif(): distribuição uniforme
rpois(): distribuição de Poisson

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

primeiro vamos ver qual é a sintaxe do comando:

?runif

agora geramos duas sequências com distribuição uniforme:

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

## [1] 9.2259964 1.6970868 5.0570761 0.4970807 5.2825136

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

## [1] 5.481742 8.644593 9.982914 6.428504 2.820517

para garantir a reprodutibilidade de simulações é sempre recomendável estabelecer uma ‘semente’ (um número inteiro)

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

A função sample() permite amostrar aleatóriamente um vetor, com ou sem substituição (replacement). É muito útil em simulações de bootstrap.

sample(1:10)

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

sample(1:10, size = 5)

## [1]  8  4  9  5 10

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

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

a função aplly(x,margin,fun) permite aplicar uma função a linhas ou colunas de matrizes:

x <- matrix(rep(1:5, 5), ncol = 5)
x

##      [,1] [,2] [,3] [,4] [,5]
## [1,]    1    1    1    1    1
## [2,]    2    2    2    2    2
## [3,]    3    3    3    3    3
## [4,]    4    4    4    4    4
## [5,]    5    5    5    5    5

# media das linhas
apply(x, 1, mean)

## [1] 1 2 3 4 5

# media das colunas
apply(x, 2, mean)

## [1] 3 3 3 3 3

# media e desvio padrão das colunas
apply(x, 2, function(y) c(mean(y), sd(y)))

##          [,1]     [,2]     [,3]     [,4]     [,5]
## [1,] 3.000000 3.000000 3.000000 3.000000 3.000000
## [2,] 1.581139 1.581139 1.581139 1.581139 1.581139

Funções semelhantes: lapply, sapply, tapply, mapply

3. Visualização de distribuições bidimensionais

Vamos simular \(n=1000\) valores de \((x,y)\) com distribuições gaussianas:

n <- 1000
x <- rnorm(mean=1, sd=1,n)
y <- rnorm(mean=2, sd=3,n)
dados <- data.frame(x,y)

# visualização
library(RColorBrewer)
library(hexbin)
rf <- colorRampPalette(rev(brewer.pal(11,'Spectral')))
r <- rf(32)

h <- hexbin(dados)
plot(h, colramp=rf)

e outro tipo de plot:

library(plot3D)

##  Create cuts:
x_c <- cut(x, 20)
y_c <- cut(y, 20)

##  Calculate joint counts at cut levels:
z <- table(x_c, y_c)

##  Plot as a 3D histogram:
hist3D(z=z, border="black")

##  Plot as a 2D heatmap:
image2D(z=z, border="black")

plot de densidades:

library(KernSmooth)

## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009

smoothScatter(x,y,nrpoints=0,add=FALSE,xlab='x',ylab='y')

cores = colorRampPalette(c("blue", "orange", "red"), space = "Lab")
smoothScatter(x,y,nrpoints=0,add=FALSE,xlab='x',ylab='y', colramp = cores)

4. exercícios/simulações:

exercícios: jogando dados

Jogamos um dado: qual é a probabilidade de o resultado ser par? A resposta, claro, é 1/2, mas vamos fazer uma simulação!

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

n=1000
# função para gerar uma distribuição uniforme e discreta entre 1 e k
rdu<-function(n,k) sample(1:k,n,replace=T)
# pdf
ddu<-function(x,k) ifelse(x>=1 & x<=k & round(x)==x,1/k,0) 
# cdf:
pdu<-function(x,k) ifelse(x<1,0,ifelse(x<=k,floor(x)/k,1))

# simulação
dado <-rdu(n,6)

# função para ver se um número é par (TRUE ou FALSE)
ehpar <- function(x){x %% 2 == 0}

par = ehpar(dado)
# número de pares:
npar = table(par)['TRUE']
# fração de pares:
npar/n

##  TRUE 
## 0.481

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=1000\) simulações onde, de cada vez, tiramos 3 bolas:

n = 1000
ncerto = 0
# o comando for() cria um loop sobre o que está dentro dos {}
for (i in 1:n){
a = sample(caixa, size = 3, replace = TRUE) 
if(a[1] == 1 & a[2] == 1 & a[3] == 3) {ncerto = ncerto + 1}
if(a[1] == 1 & a[2] == 3 & a[3] == 1) {ncerto = ncerto + 1}
if(a[1] == 3 & a[2] == 1 & a[3] == 1) {ncerto = ncerto + 1}
}
# fração simulada:
ncerto/n

## [1] 0.209

repetindo:

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

## [1] 0.178

Compare o resultado dessas duas simulações com o valor teórico. As diferenças medem 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 o comando function para fazer uma simulação com \(n\) amostras:

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

exemplos:

set.seed(1234)
simula(10)

## [1] 0.3

simula(100)

## [1] 0.17

simula(1000000)

## [1] 0.214802

Para ver esse resultado graficamente:

m=7
x=c(10,50,100,500,1000,5000,10000)
y=rep(0,m)
for (j in 1:m){
  y[j]=simula(x[j])
}
plot(x,y, main = 'tirando bolas da caixa',col='blue',xlab='n',ylab='fração')
abline(h=0.2142,col='red')

Note como a variância diminui com o número de simulações.

Para ver isso melhor, vamos calcular a variância de várias simulações (N) com \(n\) fixo.

N=100
m=4
x=rep(0,m)
x=c(10,100,1000,10000)
y=rep(0,N)
z=rep(0,m)
for (j in 1:m){
  for (k in 1:N){ y[k]=simula(x[j])}
  z[j] = sd(y)
}
# plot em log para facilitar a visualização
plot(x,z, main = 'tirando bolas da caixa',col='blue',xlab='n',ylab='sd',log='xy')
# 'ajuste' do coeficiente da curva 1/sqrt(n)
aa=mean(z*sqrt(x))
y=aa/sqrt(x)
#aa=mean(log(z))+0.5*mean(log(n))
#y=exp(aa)/sqrt(x)
lines(x,y,col='red')

A curva vermelha é \(\sigma \propto n^{-1/2}\), mostrando que a variância cai com o inverso do número de amostras.

aula 2

lsj

February 26, 2020