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- simulação do Teorema do Limite Central
  • 7- a gaussiana bivariada
  • exercícios

1 A DISTRIBUIÇÃO NORMAL OU GAUSSIANA EM R

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

distribuição gaussiana (ou normal) em R (o default é media zero e \(\sigma = 1\)):

  • 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 função de distribuição de probabilidades
  • 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

esta notação é padrão em R para as distribuições de probabilidade; exemplos:

  • dnorm, dunif, dpois, dbinon, dexp, …
  • pnorm, punif, ppois, pbinon, pexp, …
  • qnorm, qunif, qpois, qbinon, qexp, …
  • rnorm, runif, rpois, rbinon, rexp, …

densidade de probabilidades

exemplo: gaussiana de média 3 e desvio padrão 2:

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

# distribuição gaussiana
x = seq(-4, 10, by = .1)
y = dnorm(x,mean=3,sd=2)
plot(x,y, ylab='P(x)', main="função de distribuição de probabilidades gaussiana",type="l",col='blue',lwd=3)
abline(v=3)

a função dnorm() dá o valor da função de distribuição de probabilidades gaussiana em um dado ponto \(x\)

distribuição cumulativa

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

ela dá a área da função de distribuição de probabilidades de \(-\infty\) até um ponto \(x\)

# gaussiana de média 1 e desvio padrão 0.5
x = seq(-4,5,0.1)
y = pnorm(x,mean=1,sd=0.5)
plot(x,y,ylab='C(x)',main="distribuição cumulativa da gaussiana",type="l",col='gold',lwd=3)

#exemplo: a que fração da área corresponde x=-1
pnorm(-1,mean=1,sd=0.5)
## [1] 3.167124e-05
#e x=+1
pnorm(1,mean=1,sd=0.5)
## [1] 0.5
# diferença entre essas áreas (com sd=0.5, entre -sd e +sd):
pnorm(1,mean=1,sd=0.5)-pnorm(-1,mean=1,sd=0.5)
## [1] 0.4999683

vocês podem verificar que, para uma gaussiana com média 0 e sd=1,

  • a área dentro de \(\pm 1 \sigma\) é 0,68
  • a área dentro de \(\pm 2 \sigma\) é 0,95
  • a área dentro de \(\pm 3 \sigma\) é 0,997

exemplo: uma detecção é considerada real (em Física) se ela ultrapassa “5 \(\sigma\)”. Considerando uma distribuição gaussiana de média 0 e desvio padrão 1, qual é a fração da área entre -5\(\sigma\) e +5\(\sigma\)?

# o default para pnorm() é média 0 e sd = 1
# nesse caso +/-5*sigma = +/-5:
pnorm(5)-pnorm(-5)
## [1] 0.9999994

e entre 0 e 5\(\sigma\)?

pnorm(5)-pnorm(0)
## [1] 0.4999997

quantil ou percentil da distribuição cumulativa

qnorm() é o inverso de pnorm() e dá o valor de x correspondente ao percentil p da distribuição normal

p = seq(0, 1, by = .01)
q = qnorm(p)   #o default é mean=0,sd=1 

plot(p,q,main="qnorm",type="l",col='salmon',lwd=3)

# exemplo: valor de x correspondente ao percentil 0.8:
qnorm(0.8)
## [1] 0.8416212

geração de números aleatórios

rnorm() gera uma amostra com números distribuídos como uma gaussiana

# 1000 números aleatórios gerados a partir de uma distribuição gaussiana de média 10 e desvio padrão 5
# população: gaussiana de média 10 e desvio padrão 5

# dados simulados: uma amostra
x=rnorm(1000,mean=10,sd=5)

# o vetor x resultante tem 1000 elementos
# alguns valores:
head(x)
## [1]  7.197622  8.849113 17.793542 10.352542 10.646439 18.575325
# visualização
hist(x,col='salmon',breaks=20)

# media da amostra
mean(x)
## [1] 10.08064
# desvio padrão da amostra
sd(x)
## [1] 4.958475
# compare a media e o desvio padrão da amostra com o da população

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

que fração das medidas espera-se que esteja 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.375
# visualização:
hist(x,main='simulação',col='skyblue2')

# * repita o exercício com nsim = 1000000

um exemplo de modelagem gaussiana:

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

Supondo que a probabilidade da massa x (em unidades de milhões de massas solares) ser distribuída como uma gaussiana P(x)\(\sim\)N(média = 4.3, sigma= 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",col='red',lwd=3)
abline(v=5,col='blue')

probabilidade de P(x>5):

1-pnorm(5,mean=4.3,sd=0.3)
## [1] 0.009815329
# que é menor que 1%
  1. A que massa do BN corresponde o percentil de 95%?
    (a probabilidade da massa ser menor que este valor é 95%)
qnorm(0.95,mean=4.3,sd=0.3)
## [1] 4.793456
# em unidades de milhões de massas solares

A rigor, como a massa é positiva, deveríamos considerar apenas massas positivas. Verifique como isso modifica os resultados acima.

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: N=30, p=0.25

par(mfrow=c(2,2)) # panel de 2 x 2 figuras
# distribuição binomial
N=30   #número de tentativas
p = 0.25
# distribuição de probabilidades
n = seq(0,30,1) # número de sucessos - numeros inteiros positivos
y = dbinom(n,N,prob=p) 
plot(n,y,main="distribuição binomial (N,p)=(30,0.25)",type="h",lwd=2,col='green2')
# distribuição cumulativa
x=0:30
y = pbinom(x,N,prob=p) 
plot(x,y,main="distribuição cumulativa",type="h",col='green2',lwd=3)
# quantis
pp = seq(0, 1, by = .01)
q = qbinom(pp,N,prob=p)   #o default é mean=0,sd=1 
plot(pp,q,main="quantis",type="l",col='salmon',lwd=3)
# amostragem aleatória: 1000 amostras
a = rbinom(1000,N,p)
hist(a,col='salmon',xlab='n',ylab='frequência')

ilustração da distribuição binomial para vários valores de n, N e p:

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 15 objetos, qual é a probabilidade de um quasar dessa amostra ser RL?

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

dbinom(1, 15, 0.10)
## [1] 0.3431519
  1. Nessa mesma amostra, qual é a probabilidade de não mais que 2 objetos serem RL?
dbinom(0, 15, 0.10) + dbinom(1, 15, 0.10) + dbinom(2, 15, 0.10)
## [1] 0.8159389
# ou, usando a função cumulativa:
pbinom(2, 15, 0.10) 
## [1] 0.8159389

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 do número 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 dessa distribuição.

Assim, para um percentil de, digamos, 90%, com \(\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} \]

# panel de 2 x 2 figuras
par(mfrow=c(2,2))

# sequência de inteiros de 0 a 20 separados por 1
n = seq(0,20,1) # numeros inteiros positivos

# média do processo poissoniano
lambda=4

# distribuição de probabilidades
y = dpois(n,lambda) 
plot(n,y,main="dist. poissoniana (lambda=4)",type="h",lwd=2,col='blue3')

# distribuição cumulativa
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)
# quantis
y = qpois(x,lambda) 
plot(x,y,main="vetor de quantis",type="l",lwd=2,col='blue3')

# amostra poissoniana
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 qcauchy 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, \(\mu\) a média da população e \(\sigma_s\) é o desvio padrão da amostra

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", "green","red"), lty=2, cex=1)

6 Simulações do Teorema do Limite Central

6.1 simulação de distribuições exponenciais

  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] 1.210249
sd(x)
## [1] 1.154727

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
}                      
# histrogramas 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.9998899
sd(x)
## [1] 1.001613

6.2 Média e desvio padrão de uma amostra

Vamos inicialmente considerar algumas propriedades da variância: dadas duas variáveis aleatórias \(X\) e \(Y\), temos que \[V(X+Y)=V(X)+V(Y)+ Cov(X,Y)\] onde \[Cov(X,Y) = E[(X-\bar X)(Y-\bar Y)]\] é a covariância entre \(X\) e \(Y\). Se \(X\) e \(Y\) são independentes, \(Cov(X,Y)=0\).

Também: \[V(aX)=a^2 V(x)\]

Consideremos uma amostra \(\{x_1, x_2,...,x_n\}\), onde cada \(x_i\) é uma amostra independente de uma mesma distribuição de probabilidades \(P(x)\), de largura \(\sigma\), isto é, \[ \sigma^2 = V(x)\]

Se \[ \bar x = {1 \over n} \sum_i x_i, \] então, a variância da média será \[ V(\bar x) = V({1 \over n} \sum_i x_i) = {1 \over n^2}\sum_i V(x_i) = {1 \over n^2}\sum_i \sigma^2 = {n \sigma^2 \over n^2} = {\sigma^2 \over n}\] e seu desvio padrão: \[\sigma (\bar x) = {\sigma \over \sqrt{n}},\] isto é, o desvio padrão da média de uma amostra de variáveis independentes varia com o inverso da raiz quadrada do número de membros da amostra.

Vamos ver isso simulando amostras extraídas de uma distribuição uniforme:

# reprodutibilidade
set.seed(1234)
# vamos usar a biblioteca ggplot2 para fazer gráficos
library(ggplot2)

# vamos simular 1 milhão de números aleatórios entre 0 e 1
x = runif(1e6, 0, 1)

# tamanho da amostra: entre 5 e 200
n = 5:200

# função para calcular o desvio padrão da média para uma amostra de tamanho n
# a função sample() obtém n amostra do vetor x com substituição
calc_sigma_m = function(n, x) {
  sd(sample(x, n, replace = TRUE))/sqrt(n)
}

# vamos criar uma 'data frame' para armazenar os resultados para plot
df = data.frame(n,sigma_m = sapply(n, calc_sigma_m, x))

dim(df)
## [1] 196   2
head(df)
##    n    sigma_m
## 1  5 0.16523172
## 2  6 0.06937018
## 3  7 0.11902761
## 4  8 0.08217214
## 5  9 0.09912626
## 6 10 0.09558816
# figura
ggplot(df, aes(n, sigma_m)) +
  geom_line()

7. 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 2.0.0, Released 2022-12-04
## This package is based upon work supported by the National Science Foundation under Grant No. SES-0518772 and the Chan Zuckerberg Initiative: Essential Open Source Software for Science (Grant No. 2020-255193).
# 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 igual à matriz unidade em 2 dimensões

# podemos definir uma função para isso:
f <- function(x,y) dmvnorm(cbind(x,y), mu = c(0,0), sigma = diag(2))

# vejam a cara da matriz de covarância:
diag(2)
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1
# vamos criar uma grade em x e y com os valores da função em cada ponto:
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.452721 -1.241709 -0.18939477
## [2,]  1.457649 -1.825032  0.07419385
## [3,] -0.893978  1.554406 -0.07609819
# 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)
# alguns pontos simulados:
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
# coordenadas dos pontos simulados
X = gbv[,1]
Y = gbv[,2]

# visualização
# função para desenhar elipses para a distribuição normal bivariada
# 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)

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

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 a) se observar menos que 15 objetos e b) 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 gaussiana bivariada de média \(\mathbf \mu =\{0,0\}\) e variâncias \(\sigma_x = 1\) e \(\sigma_y = 2\). Adote para o coeficiente de correlação \(\rho = 0.7+ g\), onde \(g\) é o valor resultante de: \(set.seed(SeuNumeroUsp); g = 0.3*runif(1)\). Determine a distribuição de probabilidade condicional \(P(y|x=1)\).