aula 2
lsj
December 27, 2018
aula de hoje: Introdução ao R (II): probabilidades
distribuição gaussiana
simulações
visualização de distribuições bidimensionais
valores esperados
exercícios/simulações
1. DISTRIBUIÇÃO GAUSSIANA
distribuição gaussiana de média \(\mu\) e variância \(\sigma^2\): \[ 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
hist(rnorm(100,mean=10,sd=5))
Exemplo:
- 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)?
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.7934562. 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 sguidas de 5 números com distribuição uniforme entre 0 e 10:
primeiro vamos ver qual é a sintaxe do comando:
?runifagora geramos duas sequências com distribuição uniforme:
runif(5,min=0,max=10)## [1] 0.4275082 3.2663596 3.3581583 5.8252975 3.9979338runif(5,min=0,max=10)## [1] 3.7988215 7.6822267 8.2605206 0.4940009 4.3562918para garantir a reprodutibilidade de simulações é 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.609154set.seed(1234)
runif(5,min=0,max=10)## [1] 1.137034 6.222994 6.092747 6.233794 8.609154A 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] 7 1 2 5 4 6 3 9 8 10sample(1:10, size = 5)## [1] 9 3 10 2 7sample(1:10, size = 15, replace = TRUE)## [1] 4 4 2 1 3 9 6 10 9 1 5 3 4 6 2a 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 5apply(x, 1, mean)## [1] 1 2 3 4 5apply(x, 2, mean)## [1] 3 3 3 3 3apply(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.581139Funçõ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)
df <- data.frame(x,y)
# visualização
library(RColorBrewer)
library(hexbin)
rf <- colorRampPalette(rev(brewer.pal(11,'Spectral')))
r <- rf(32)
h <- hexbin(df)
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")
4. Valores esperados
valores esperados de \(x\) e \(y\) (os valores simulados foram 1 e 2):
mean(df$x)## [1] 1.002755mean(df$y)## [1] 2.065783# ou usando apply
apply(df,2,mean)## x y
## 1.002755 2.065783valores esperados dos desvios padrão (os valores simulados foram 1 e 3:
sd(df$x)## [1] 0.9957839sd(df$y)## [1] 3.040505# ou usando apply
apply(df,2,sd)## x y
## 0.9957839 3.0405046outro tipo de plot:
library(KernSmooth)## KernSmooth 2.23 loaded
## Copyright M. P. Wand 1997-2009smoothScatter(x,y,nrpoints=0,add=FALSE,xlab='x',ylab='y',main = 'gaussiana bidimensional')
cores = colorRampPalette(c("blue", "orange", "red"), space = "Lab")
smoothScatter(x,y,nrpoints=0,add=FALSE,xlab='x',ylab='y',main = 'gaussiana bidimensional', colramp = cores)
5. 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.508exercí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.2repetindo:
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.177Compare 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.3simula(100)## [1] 0.29simula(1000000)## [1] 0.214924Para 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.