aula de hoje: distribuições de probabilidades

  • 1 distribuição gaussiana
  • 2 distribuição binomial
  • 3 distribuição de Poisson
  • 4 funções de distribuição em R
  • 5 combinação de distribuições de probabilidades
  • 6 simulação do Teorema do Limite Central por Monte Carlo

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")
# distribuição cumulativa
plot(x,pnorm(x,mean=0,sd=1),main="distribuição cumulativa",type="l")
# percentil da distribuição cumulativa
p = seq(0, 1, by = .01)
x = qnorm(p)
plot(p,x,main="função quantil",type="l")
# geração de números aleatórios
hist(rnorm(100,mean=10,sd=5),main='simulação')

2 DISTRIBUIÇÃO BINOMIAL

probabilidade de n sucessos em N tentativas independentes onde a probabilidade de êxito em cada tentativa é a mesma e vale \(\theta\): \[ P(n) = \binom Nn \theta^n (1-\theta)^{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")
y = dbinom(n,N,prob=0.10) 
plot(n,y,main="distribuição binomial (N,p)=(30,0.1)",type="h")

# 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")
N=30
y = dbinom(n,N,prob=0.5) 
plot(n,y,main="distribuição binomial (N,p)=(30,0.5)",type="h")

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 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")
y = ppois(n,lambda) 
plot(n,y,main="distr. cumulativa",type="s")
# 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")
x=rpois(100,lambda)
hist(x,main="distr. aleatória com lambda=4")

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 pelo menos 20 quasares por unidade de área?

a probabilidade de termos 19 quasares ou menos é

ppois(19, 15) #lower tail
## [1] 0.8752188

a probabilidade de 20 ou mais quasares por unidade de área é

ppois(19, 15, lower=FALSE) #upper tail
## [1] 0.1247812

que é equivalente a

1-ppois(19, 15)
## [1] 0.1247812

4 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))
lines(x,y2,type='l',col='blue')
lines(x,y3,type='l',col='red')
legend(2, 0.5, legend=c("dof=2", "dof=4", "dof=6"),
       col=c("black", "blue","red"), lty=1, cex=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)

5 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. 
## -1767371.0       -0.8        0.1    -1765.9        1.1      704.5

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

6 Simulação do Teorema do Limite Central por Monte Carlo

  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")

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

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