aula de hoje: distribuições de probabilidades

  • 1 distribuição gaussiana
  • 2 distribuição binomial
  • 3 distribuição exponencial
  • 4 distribuição de Poisson
  • 5 funções de distribuição em R
  • 6 combinação de distribuições de probabilidades
  • 7 simulação do Teorema do Limite Central
  • 8 a gaussiana bivariada
  • 9 exercícios

1 A DISTRIBUIÇÃO GAUSSIANA EM R

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- quartis: inverso da função de distribuição cumulativa

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

plots:

# reprodutibilidade
set.seed(1234)

# densidade de probabilidades, distribuição cumulativa e percentil da distribuição cumulativa
par(mfrow=c(1,4)) # 4 figuras lado a lado
x = seq(-4, 4, by = .1) # sequência de pontos para x
# densidade de probabilidades
plot(x,dnorm(x,mean=0,sd=1),main="pdf gaussiana", type="l", lwd=2,col='darkgoldenrod1')
# distribuição cumulativa
plot(x,pnorm(x,mean=0,sd=1),main="distribuição cumulativa",type="l",lwd=2,col='darkgoldenrod1')
# percentil da distribuição cumulativa
p = seq(0, 1, by = .01)
x = qnorm(p)
plot(p,x,main="função quantil",type="l",lwd=2,col='darkgoldenrod1')
# geração de números aleatórios
hist(rnorm(100,mean=10,sd=5),main='simulação',col='salmon')

Exemplo: suponha que temos medidas de magnitudes de uma estrela em uma certa banda fotométrica e que pode ser modelada como uma gaussiana de média \(\mu = 20\) e desvio padrão \(\sigma = 1\)

  • que fração das medidas espera-se que estejam dentro de 0.5 magnitude da média?
# parâmetros da gaussiana: média e desvio padrão
mu = 20
sigma = 1
pnorm(mu+0.5,mean=mu,sd=sigma)-pnorm(mu-0.5,mean=mu,sd=sigma)
## [1] 0.3829249
# vamos fazer uma simulação:
nsim = 1000
x = rnorm(nsim,mean=mu,sd=sigma)
# fração dos pontos no intervalo desejado:
length(x[x > mu-0.5 & x < mu+0.5])/nsim
## [1] 0.4
# * repita o exercício com nsim = 1000000

2 DISTRIBUIÇÃO BINOMIAL

probabilidade de n sucessos em N tentativas independentes onde a probabilidade de êxito em cada tentativa é a mesma e vale \(p\): \[ P(n) = \binom Nn p^n (1-p)^{N-n}\] distribuição binomial em R:

dbinom(x, size, prob, log = FALSE)

pbinom(q, size, prob, lower.tail = TRUE, log.p = FALSE)

qbinom(p, size, prob, lower.tail = TRUE, log.p = FALSE)

rbinom(n, size, prob)

ilustração da distribuição binomial:

par(mfrow=c(2,2)) # panel de 2 x 2 figuras
# distribuição binomial com p=0.1
N=15   #número de tentativas
n = seq(0,20,1) # número de sucessos - numeros inteiros positivos
y = dbinom(n,N,prob=0.10) 
N=30
plot(n,y,main="distribuição binomial (N,p)=(15,0.1)",type="h",lwd=2,col='green2')
y = dbinom(n,N,prob=0.10) 
plot(n,y,main="distribuição binomial (N,p)=(30,0.1)",type="h",lwd=2,col='green2')

# distribuição binomial com p=0.5
N=15
y = dbinom(n,N,prob=0.5) 
plot(n,y,main="distribuição binomial (N,p)=(15,0.5)",type="h",lwd=2,col='green2')
N=30
y = dbinom(n,N,prob=0.5) 
plot(n,y,main="distribuição binomial (N,p)=(30,0.5)",type="h",lwd=2,col='green2')

propriedades da distribuição binomial:

  • média = Np;
  • variância = Np(1-p);

Compare com a figura acima.

Exemplos:

  1. Em uma amostra de quasares, espera-se que 10% sejam "radio-loud" (RL). Numa amostra de 25 objetos, qual é a probabilidade de um quasar dessa amostra ser RL?

Isso pode ser rescrito como dbinom(1, 25, 0.10):

dbinom(1, 25, 0.10)
## [1] 0.1994161
  1. Nessa mesma amostra, qual é a probabilidade de não mais que 2 objetos serem RL?
pbinom(2, 25, 0.10) #função cumulativa
## [1] 0.5370941

3 DISTRIBUIÇÃO EXPONENCIAL

Suponha que o núcleo de um dado elemento químico radioativo tenha uma meia vida distribuida exponencialmente com taxa \(\lambda\). Seja \(n(t)\) o número esperado de núcleos no instante \(t\), com \(n(t=0)=n_0\). Qual é o desvio padrão esperado em medidas de \(n(t)\)?

modelo probabilístico: a probabilidade de um núcleo decair no instante \(t\) é \[ P(t|\lambda) = \left \{ \begin{matrix} \lambda e^{-\lambda t}, & \mbox{se }t \ge 0 \\ 0, & \mbox{se } t<0 \end{matrix} \right. \] média e variância: \[ \bar t = \frac{1}{\lambda} ~~~~~~\sigma^2 = \frac{1}{\lambda^2}\]

valor esperado de núcleos no instante \(t\): \(n(t) = n_0 e^{-\lambda t}\)

desvio padrão: \(n(t)/\lambda\)

vamos fazer um estudo numérico para verificar:

# reprodutibilidade
set.seed(321)

# meia vida (em unidades arbitrárias)
meiavida = 2

# taxa: decaimentos por unidade de tempo
lambda = 1/meiavida
lambda
## [1] 0.5
# número inicial de átomos
n0 = 1000

# tempos de decaimento de cada átomo
td = rexp(n0,rate=lambda)
summary(td)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
##  0.001222  0.574221  1.360260  2.005711  2.973120 16.655094
hist(td,breaks=20, col='salmon', main='tempos de decaimento')

# número de átomos sobreviventes em um tempo t
# vamos usar a função "empirical cumulative distribution function":
nt = ecdf(td)
plot(nt,lwd=2,col='cadetblue')

# probabilidade de que um núcleo sobreviva sem decair por pelo menos 10 unidades de tempo:
nn = sort(td)
p10 = length(nn[nn > 10])/n0
p10
## [1] 0.004
# valor teórico: exp(-lambda*10)=exp(-5)=
exp(-5)
## [1] 0.006737947
# vamos simular 1000 vezes esse processo
nsim = 1000
p10 = NULL
for(i in 1:nsim){
# tempos de decaimento
td = rexp(n0,rate=lambda)

# número de átomos sobreviventes em um tempo t
# vamos usar a função "empirical cumulative distribution function":
nt = ecdf(td)
if(i == 1)plot(nt)
lines(nt,col='red')

# probabilidade de que um núcleo sobreviva sem decair por pelo menos 10 unidades de tempo:
nn = sort(td)
p10[i] = length(nn[nn > 10])/n0
}
x=seq(0,15,0.01)
y=1-exp(-lambda*x)
lines(x,y,lwd=2)

hist(p10,breaks=20, col='salmon', main='')
abline(v=exp(-5),col='green',lwd=2)
abline(v=mean(p10),col='blue',lwd=2)
legend('topright', legend=c("simulação", "valor teórico"),
       col=c("blue","green"), lty=1, cex=1,,box.lty=0,lwd=2)

Para a distribuição exponencial é fácil determinar percentis.

A distribuição cumulativa da distribuição exponencial é \[ F(x) = \int_0^x \lambda e^{-\lambda t} dt = 1-e^{-\lambda x}\]

O percentil \(p\) corresponde ao valor de x tal que \(F(x)=p\) e, portanto, \[1-e^{-\lambda x} = p\] ou \[x = -{1 \over \lambda}log(1-p). \] que é a função quantil.

Assim, para um percentil de, digamos, 90%, e \(\lambda=2\), temos que

x = -1/2*log(1-0.9)
x
## [1] 1.151293

4 DISTRIBUIÇÃO DE POISSON

probabilidade de \(n\) eventos ocorrerem num certo intervalo de tempo, ou em certa região do espaço, se estes eventos ocorrem independentemente e com uma média fixa \(\lambda\): \[ P(n) = {\lambda^n \over n!} e^{-\lambda} \]

par(mfrow=c(2,2)) # panel de 2 x 2 figuras
# sequência de inteiros de 0 a 20 separados por 1
n = seq(0,20,1) # numeros inteiros positivos
lambda=4
y = dpois(n,lambda) 
plot(n,y,main="dist. poissoniana (lambda=4)",type="h",lwd=2,col='blue3')
y = ppois(n,lambda) 
plot(n,y,main="distr. cumulativa",type="s",lwd=2,col='blue3')
# 100 números uniformemente distribuídos entre 0 e 1
x = seq(0,1,0.001)
y = qpois(x,lambda) 
plot(x,y,main="vetor de quantis",type="l",lwd=2,col='blue3')
x=rpois(100,lambda)
hist(x,main="distr. aleatória com lambda=4",lwd=2,col='skyblue')

Exemplos:

  1. Se em um certo intervalo de redshift a densidade de quasares é de 15 por unidade de área, qual a probabilidade de um certo campo ter a) menos de 20 quasares e b) pelo menos 20 quasares por unidade de área?
  1. a probabilidade de termos 19 quasares ou menos é
ppois(19, 15) #lower tail
## [1] 0.8752188
  1. a probabilidade de 20 ou mais quasares por unidade de área é
ppois(20, 15, lower=FALSE) #upper tail
## [1] 0.08297091

que é equivalente a

1-ppois(20, 15)
## [1] 0.08297091
  1. Dois contadores Geiger registram a chegada de raios cósmicos, um com uma taxa de 3 eventos por segundo e o outro com 4 eventos por segundo. Suponha que o processo seja poissoniano. Num certo segundo, os detectores tiveram 8 contagens. Qual é a probabilidade de que cada contador detectou 4 eventos?

Como os detectores e as medidas são independentes, sendo \(P(n|\mu)\) a probabilidade de se obter \(n\) contagens num processo poissoniano de média \(\mu\) \[ P(n=4|\mu=3).P(n=4|\mu=4) \simeq 0.0328 \]

em R:

dpois(4,lambda=3)*dpois(4,lambda=4)
## [1] 0.03282775

Vamos verificar fazendo uma simulação:

# reprodutibilidade
set.seed(456)

# número de simulações
nsim=100000

# simulação das contagens do contador 1
n1 = rpois(nsim,lambda=3)
# alguns valores
n1[1:10]
##  [1] 1 2 4 5 4 2 1 2 2 2
# simulação das contagens do contador 2
n2 = rpois(nsim,lambda=4)
n2[1:10]
##  [1] 8 3 4 4 2 4 3 7 4 4
# vamos selecionar os casos em que n1 e n2 são iguais a 4
sel = rep(0,nsim)
sel[n1 == 4 & n2 == 4] = 1
sum(sel)
## [1] 3194
# fração das simulações com n1 e n2 iguais a 4
sum(sel)/nsim
## [1] 0.03194
# compare esse valor com o esperado teoricamente

Um resultado importante: se \(X_1\) e \(X_2\) são variáveis com distribuição poissoaniana com médias \(\mu_1\) e \(\mu_2\), respectivamente, a distribuição de \(X_1+X_2\) será poissoniana de média \(\mu_1+\mu_2\).

No exemplo acima, \(\mu_1=3\) e \(\mu_2=4\) e \(\mu_1+\mu_2=7\)

Compare com o valor obtido na simulação:

mean(n1+n2)
## [1] 7.01165

5 FUNÇÕES DE DISTRIBUIÇÃO EM R

nome p q d r
Beta pbeta qbeta dbeta rbeta
Binomial pbinom qbinom dbinom rbinom
Cauchy pcauchy cauchy dcauchy rcauchy
Chi-Square pchisq qchisq dchisq rchisq
Exponential pexp qexp dexp rexp
F pf qf df rf
Gamma pgamma qgamma dgamma rgamma
Geometric pgeom qgeom dgeom rgeom
Hypergeometric phyper qhyper dhyper rhyper
Logistic plogis qlogis dlogis rlogis
Log Normal plnorm qlnorm dlnorm rlnorm
Negative Binomial pnbinom qnbinom dnbinom rnbinom
Normal pnorm qnorm dnorm rnorm
Poisson ppois qpois dpois rpois
Student t pt qt dt rt
Studentized Range ptukey qtukey dtukey rtukey
Uniform punif qunif dunif runif
Weibull pweibull qweibull dweibull rweibull
Wilcoxon Rank Sum Statistic pwilcox qwilcox dwilcox rwilcox
Wilcoxon Signed Rank Statistic psignrank qsignrank dsignrank rsignrank

exemplo: distribuição do \(\chi^2\)

Suponha que temos dados \(\{x_i\}\), \(i=1,n\), amostrados com \(N(\mu,\sigma)\).

Sendo \[z_i = \frac{x_i-\mu}{\sigma}\] (os resíduos normalizados) e \[Q = \sum_{i=1}^n z_i^2,\] (a soma do quadrado desses resíduos) então \(Q\) segue uma distribuição de \(\chi^2\) com \(k=n\) graus de liberdade: \[ P(Q|k) = \chi^2(Q|k) = \frac{1}{2^{k/2}\Gamma(k/2)} Q^{\frac{k}{2}-1} e^{-\frac{Q}{2}} \]

Exemplos:

x=seq(0,6,0.01)
y1=dchisq(x,2)
y2=dchisq(x,4)
y3=dchisq(x,6)
plot(x,y1,type='l',ylab=expression(chi^2),lwd=2)
lines(x,y2,type='l',col='blue',lwd=2)
lines(x,y3,type='l',col='red',lwd=2)
legend(2, 0.5, legend=c("dof=2", "dof=4", "dof=6"),
       col=c("black", "blue","red"), lty=1, cex=2,lwd=2)

exemplo: distribuição \(t\) de Student

Se temos dados \(\{x_i\}\), \(i=1,n\), amostrados com \(N(\mu,\sigma)\), então a variável \[ t = \frac{\bar x - \mu}{\sigma_s/\sqrt{N}} \] se distribui com uma pdf \(t\) com \(k=n-1\) graus de liberdade.

\(\bar x\) é a média da amostra e \(\mu\) a média da população

Note que \(t\) é baseado nas estimativas \(\bar x\) e \(\sigma_s\), enquanto que o \(\chi^2\) é baseado nos valores verdadeiros \(\mu\) e \(\sigma\)

Exemplos:

x=seq(-5,5,0.1)
y=dnorm(x)
y1=dt(x,1)
y2=dt(x,2)
y3=dt(x,10)
plot(x,y,type='l',ylab='distribuição t',main= 'gaussiana em azul',col = 'blue')
lines(x,y1,type='l',col='black',lty=2)
lines(x,y2,type='l',col='green',lty=2)
lines(x,y3,type='l',col='red',lty=2)
legend(-5, 0.4, legend=c("dof=1", "dof=2", "dof=10"),
       col=c("black", "blue","red"), lty=2, cex=1)

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

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) = \frac{1}{2 \pi}~~~~ (0 < \phi < 2\pi) \] e que \[ y = \sin (\phi). \] 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}} \]

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

n=100000L  #esse L significa que n é inteiro
fi=runif(n, min = 0, max = 2*pi)
y=sin(fi)
h=hist(y)

hmin=min(h$counts) #para normalizar P(y)
y=seq(-0.99,0.99,0.01)
yc=hmin/sqrt(1-y^2)
lines(y,yc,col='red',lwd=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 distruição de \(z\) fazendo uma simulação:

par(mfrow=c(1,2))
m=1000
x=rnorm(m)
y=rnorm(m)
z=x/y
hist(z,main="distribuicao de z=x/y")
zz = z[z > -5 & z < 5]
hist(zz,main="distribuicao de z=x/y",xlim=c(-5,5),breaks=100)

# note os outliers!
summary(z)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## -960.1552   -0.8250    0.0038    0.4894    0.9963 1742.5081

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

par(mfrow=c(2,2)) # panel de 2 x 2 figuras
# sequência de inteiros de 0 a 20 separados por 1
x = seq(-5,5,0.1) # numeros inteiros positivos
y = dcauchy(x) 
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)
y = pcauchy(x) 
plot(x,pnorm(x),main="distr. cumulativa",type="l",col='red')
lines(x,y,col='blue')
# 100 números uniformemente distribuídos entre 0 e 1
x = seq(0,1,0.01)
y = qcauchy(x) 
plot(x,y,main="vetor de quantiles",type="l",col='blue')
lines(x,qnorm(x),col='red')
x=rcauchy(100)
hist(x,main="distr. aleatória ",col='blue')

7 Simulação do Teorema do Limite Central

  1. Vamos simular n conjuntos de m pontos gerados com uma distribuição exponencial
  1. Para cada simulação calculamos a média dos pontos
  1. Plotamos as médias das distribuições para verificar se a distribuição das médias se assemelha a uma gaussiana

Inicialmente, vamos gerar m=10 pontos distribuídos exponencialmente e calcular a média

m=10
x=rexp(m)
y=exp(-x)
plot(x,y,main="uma simulação com 10 pontos",pch=20,col='red',cex=2)

# média e desvio padrão desses pontos:
mean(x)
## [1] 0.8522886
sd(x)
## [1] 1.279011

Agora um exemplo de programa para esta simulação:

par(mfcol=c(2,2))                   # definimos a forma de apresentação das figuras
m=10                                # número de amostras de cada exponencial
for (n in c(10,100,1000,100000))   # loop (externo) da simulação: 4 conjuntos de n
{
xm =rep(0,n)                        # crio um vetor de dimensão n com zeros
for (i in 1:n)                      #loop (interno) da simulacao sobre m
{
x=rexp(m)                           # gero m pontos com distribuição exponencial        
xm[i]=mean(x)                       # guardo a média dos m pontos no vetor xm
}                      # estatísticas de xm
hist(xm,main=paste("Nsim = ",n),col='salmon')    # histogramas de xm
}

Valores esperados:

  • média: \(\bar x = \int_0^{\infty} x P(x) dx = \int_{0}^{\infty} x e^{-x} dx = 1\)
  • desvio padrão: \(\sigma = \Bigg( \int_{0}^{\infty} (x - \bar x)^2 e^{-x} dx \Bigg)^{1/2} = 1\)

Vamos verificar:

m=1000000
x=rexp(m)
y=exp(-x)
# média e desvio padrão desses pontos:
mean(x)
## [1] 0.9994695
sd(x)
## [1] 1.001051

Exemplo: random walk em uma dimensão

Já vimos que a média de uma variável se distribui como \[\bar x \sim N(\mu, \sigma/\sqrt{N})\]

  • Seja \(S = x_1+x_2+...+x_n\) a soma de \(n\) variáveis aleatórias independentes extraídas da mesma distribuição \(P(x)\), que tem média \(\mu\) e desvio padrão \(\sigma\).
  • O teorema do limite central diz que, para \(n\) grande, a distribuição de \(S\) é aproximadamente gaussiana, de média \(\bar x = S/n\) e desvio padrão da média \(\hat \sigma = \sigma / \sqrt{n}\).
  • O random walk é usado para modelar a difusão ou os movimentos aleatórios de partículas
  • Considere uma partícula movendo-se ao longo do eixo \(x\) com deslocamentos, a partir de um dado ponto \(x_i\), de valor -1, 0 ou +1, todos com probabilidade 1/3.
  • A posição \(S\) de uma partícula no tempo \(n\) é a soma dos deslocamentos \(x_1,...,x_n\).
  • Problema: qual é a probabilidade de, após 10000 passos, a partícula estar a mais de 100 unidades do ponto inicial (\(|S| > 100\))?
  • Seja \(x\) um único passo: \[\mu = E(x) = {-1 \over 3}+{0 \over 3}+{1 \over 3} = 0 \] \[ V(x)= {(-1)^2 \over 3}+{0^2 \over 3}+{1^2 \over 3} = {2 \over 3} \] \[ \sigma = \sqrt{2/3} \simeq 0.8165\] assim, \[ E(S) = 10000 \mu = 0\] \[ \sigma(S) = \sqrt{10000}.\sigma \simeq 81.65\] e a probabilidade de se ter \(x>100\) em uma gaussiana de média 0 e desvio padrão 81.65 é:
1-pnorm(100,mean=0,sd=81.65)
## [1] 0.1103366
  • A probabilidade de se ter \(|x|>100\), é duas vezes esse valor, \(\sim 22\)%.
  • Vamos verificar isso com simulação: vamos simular 1000 vezes o problema
# número de simulacões
nsim = 1000
# número de passos em cada simulação
npassos = 10000
# inicialização de variáveis
n100 = 0
X=NULL
S=NULL
# simulações
for(k in 1:nsim){
# passos
DX = sample(c(-1,0,1), size = npassos, replace = TRUE)
# começa-se da origem
X[1]=0
for(i in 2:npassos){X[i]=X[i-1]+DX[i]}
S[k]=X[npassos]
if(abs(S[k]) >= 100)n100=n100+1
if(k == 1)hist(X,main='exemplo de uma simulação',breaks=20,col='salmon')
}

hist(S,main='',breaks=20,col='salmon')

# probabilidade:
n100/nsim
## [1] 0.236

8. A distribuição gaussiana bivariada

\[ P(x,y|\mu_x, \mu_y, \sigma_x, \sigma_y, \sigma_{xy}) = \frac{1}{2 \pi \sigma_x \sigma_y \sqrt{1-\rho^2}} \exp \Bigg(-\frac{Z^2}{2(1-\rho^2)} \Bigg) \] onde \[ Z^2 = \frac{(x-\mu_x)^2}{\sigma_x^2} + \frac{(y-\mu_y)^2}{\sigma_y^2} -2\rho \frac{(x-\mu_x)(y-\mu_y)}{\sigma_x \sigma_y} \] coeficiente de correlação: \[ \rho = \frac{\sigma_{xy}}{\sigma_x \sigma_y} \]

Em R temos vários pacotes para simular a gaussiana bivariada. Aqui vamos usar o pacote mixtools, onde as funções dmvnorm() e rmvnorm() são semelhantes a dnorm() e rnorm():

#library("mvtnorm")
library('mixtools')
## mixtools package, version 1.1.0, Released 2017-03-10
## This package is based upon work supported by the National Science Foundation under Grant No. SES-0518772.
# vamos gerar 30 pontos em x e em y entre -3 e 3
x <- y <- seq(from = -3, to = 3, length.out = 30)

# vamos plotar uma gaussiana bivariada de media mu=(0,0) e matriz de covariância a matriz unidade em 2 dimensões
#f <- function(x,y) dmvnorm(cbind(x,y), mean = c(0,0), sigma = diag(2))
f <- function(x,y) dmvnorm(cbind(x,y), mu = c(0,0), sigma = diag(2))

diag(2)
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1
# criamos uma grade em x e y
z <- outer(x, y, FUN = f)
# plot
persp(x, y, z, theta = -30, phi = 30, ticktype = "detailed",col='gold')

# o pacote funciona para gaussianas multivariadas em geral!
# ex: amostragem de uma gaussiana multivariada em 3 dimensões
# simulação  de n pontos
n=3
rmvnorm(n,mu = c(0,0,0), sigma = diag(3))
##            [,1]      [,2]       [,3]
## [1,]  1.4039554  1.311689  0.8687371
## [2,] -0.7145790 -0.066788  1.4960679
## [3,]  0.6621682 -1.182961 -0.2488582
# façam ?dmvnorm, ?outer e ?persp

Suponha que \(u\) e \(v\) tenham uma distribuição gaussiana bivariada com parâmetros \(\mu_u,\sigma_u,\mu_v,\sigma_v\). Muitas vezes é conveniente "padronizar" (standardize) as variáveis: \[ x = {u -\mu_u \over \sigma_u}~~~~~~~~~y= {v -\mu_v \over \sigma_v} \] Pode-se verificar que \[\rho = Cov(u,v) = Cov(x,y)\]

Vamos simular \(N\) amostras de uma gaussiana bivariada:

# reprodutibilidade
set.seed(123)
library(mixtools)
# número de amostras
N = 2000
# parâmetros da distribuição
rho = -0.6
mu1 = 1
s1 = 2
mu2 = 1
s2 = 8

# notação vetorial
# média
media = c(mu1,mu2)
# matriz de covariância
Sigma = matrix(c(s1^2, s1*s2*rho, s1*s2*rho, s2^2),nrow=2,ncol=2)
Sigma
##      [,1] [,2]
## [1,]  4.0 -9.6
## [2,] -9.6 64.0
# simulação
gbv = rmvnorm(N, media, Sigma)
head(gbv)
##            [,1]        [,2]
## [1,]  0.2551418 -0.27105266
## [2,]  3.6367117  0.01303956
## [3,] -0.4773147 14.48627899
## [4,]  3.0556669 -9.49969516
## [5,]  0.2495100 -1.85621522
## [6,]  2.7685660  2.64165312
X = gbv[,1]
Y = gbv[,2]

# visualização
# Function to draw ellipse for bivariate normal data
# https://blog.revolutionanalytics.com/2016/08/simulating-form-the-bivariate-normal-distribution-in-r-1.html
# a função ellipse() é do pacote mixtools 
ellipse_bvn <- function(bvn, alpha){
  Xbar <- apply(bvn,2,mean)
  S <- cov(bvn)
  ellipse(Xbar, S, alpha = alpha, col="red")
}
# alpha é o nível de confiança:
#0.5- espera-se que 50% dos pontos estarão dentro da elipse
#0.05- espera-se que 95% dos pontos estarão dentro da elipse

par(mfrow=c(1,1))
plot(X,Y, xlab="X",ylab="Y",main= "distribuição gaussiana bivariada",col='grey50')
ellipse_bvn(gbv,.5)
ellipse_bvn(gbv,.05)

# distribuições marginais
par(mfrow=c(1,2))
hist(X,xlab= 'X',main = 'distribuição de X',col='salmon')
abline(v = mean(X))
abline(v = mu1, col='blue')
hist(Y,xlab= 'Y',main = 'distribuição de Y',col='salmon')
abline(v = mean(Y))
abline(v = mu2, col='blue')

# distribuição marginal de X para 7 < Y < 10
par(mfrow=c(1,1))
hist(X[Y > 7 & Y < 10],xlab= 'X',main = '7 < Y < 10',col='salmon')

# vamos comparar as médias e desvios padrão da distribuição e da amostra:
c(mu1,mean(X))
## [1] 1.000000 1.017071
c(mu2,mean(Y))
## [1] 1.000000 1.012898
c(s1,sd(X))
## [1] 2.00000 2.00843
c(s2,sd(Y))
## [1] 8.000000 8.014775
# coeficiente de correlação entre as duas variáveis:
cor(X,Y)
## [1] -0.6209815

Pode-se simular a distribuição gaussiana bivariada a partir de distribuições gaussianas unidimensionais usando a probabilidade condicional de \(y\) para um dado \(x\):
\[ P(y|x) \sim N(\mu_{y|x},\sigma_{y|x}), \] onde \[ \mu_{y|x} = \mu_y + \rho \frac{\sigma_y}{\sigma_x} (x-\mu_x) \] \[ \sigma_{y|x} = \sigma_y\sqrt{1-\rho^2} \]

Nesse caso simula-se \(x\) como \(P(x) \sim N(\mu_x, \sigma_x)\) e, para cada \(x\), simula-se \(y\) como \(P(y|x)\):

#  função que amostra a distribuição
rgbv<-function (n, rho, mu1, s1, mu2, s2) 
{
        x <- rnorm(n, mean=mu1, sd=s1)
        y <- rnorm(n, mu2+rho*s2/s1*(x-mu1), s2*sqrt(1 - rho^2))
        cbind(x, y)
}

# simulação 
gbv<-rgbv(N,rho, mu1, s1, mu2, s2)

# visualização
par(mfrow=c(1,1))
X = gbv[,1]
Y = gbv[,2]
plot(X,Y, xlab="x",ylab="y",main= "distribuição gaussiana bivariada",col='grey60')
ellipse_bvn(gbv,.5)
ellipse_bvn(gbv,.05)

Esee tipo de amostragem, usando probabilidades condicionais, é chamado de amostragem de Gibbs.

9. Exercícios

    1. Se uma fração \(f\) das estrelas são binárias, qual é a probabilidade de que uma amostra aleatória de 5 estrelas contenha a) 0, b) 1, c) 2, d) 3, e) 4 e f) 5 binárias?
      Adote para \(f\) o valor resultante de:
      \(set.seed(SeuNumeroUsp); f = 0.5+0.2*runif(1)\)
      Represente graficamente seu resultado.
    1. Em um levantamento de quasares cobrindo 1000 graus quadrados no céu encontrou-se \(N\) quasares. Se a distribuição projetada destes objetos for modelada como uma distribuição de Poisson, qual é a probabilidade de se observar menos que 15 objetos e mais que 30 em uma certa área de 1 grau quadrado?
      Adote para \(N\) o valor resultante de: \(set.seed(SeuNumeroUsp); N = as.integer(20000+5000*runif(1))\)
      Teste seu resultado com simulações.
    1. Considere a distribuição \(P(x,y) = a(x-y)^2\), com \(0<x<1\) e \(0<y<1\)
      1. determine \(a\) para que \(P(x,y)\) seja uma distribuição de probabilidades
      1. determine \(P(x)\) e \(P(y)\)
    1. Considere a distribuição gaussiana bivariada de média \(\mathbf \mu =\{0,0\}\), variâncias \(\sigma_x = 1\), \(\sigma_y = 2\) e coeficiente de correlação \(\rho = 0.7\). Determine a distribuição de probabilidade condicional \(P(y|x=1)\).