distribuições de probabilidades
lsj
February 26, 2020
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
7 gaussiana bivariada
8 ajuste de uma gaussiana aos dados
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- 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)) # 3 figuras em uma linha
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 da gaussiana",type="l")
# percentil da distribuição cumulativa
p = seq(0, 1, by = .01)
x = qnorm(p)
plot(p,x,main="qnorm",type="l")
# geração de números aleatórios
hist(rnorm(100,mean=10,sd=5),main='simulação de uma gaussiana')
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")
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 = Npq;
Compare com a figura acima.
Exemplos:
- 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- 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.53709413 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 quantiles",type="l")
x=rpois(100,lambda)
hist(x,main="distr. aleatória com lambda=4")
Exemplos:
- 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.8752188a probabilidade de 20 ou mais quasares por unidade de área é
ppois(19, 15, lower=FALSE) #upper tail## [1] 0.1247812que é equivalente a
1-ppois(19, 15)## [1] 0.12478124 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}\] e \[Q = \sum_{i=1}^n z_i^2,\] 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.
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 (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
Vamos 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:
m=1000
x=rnorm(m)
y=rnorm(m)
z=x/y
hist(z,main="distribuicao de z=x/y")
summary(z)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -309.75458 -0.81996 0.03631 -0.52386 1.11215 190.98368Pode-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
Vamos simular n conjuntos de m pontos gerados com uma distribuição exponencial
Para cada simulação calculamos a média dos pontos
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] 1.231798sd(x)## [1] 1.131787Agora 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.001172sd(x)## [1] 1.0031857 Simulação de uma gaussiana bivariada
Para criar amostras distribuídas como uma gaussiana bivariada vamos usar o algoritmo de Gibbs, onde as amostras em \(x\) e \(y\) são obtidas de forma sucessiva e independente.
A gaussiana bivariada pode ser definida como \[ 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} \] e \(\rho\) é o coeficiente de correlação: \[ \rho = \frac{\sigma_{xy}}{\sigma_x \sigma_y} \]
Pode-se mostrar que a distribuição univariada de \(x\) para \(y\) fixo é (ver completando_o_quadrado.pdf) é \[P(x|y) \sim N(\mu,\sigma)\] onde \[ \mu = \mu_x - \frac{\rho \sigma_x (y-\mu_y)}{\sigma_y}\] e \[ \sigma = \frac{\sigma_x}{(1-\rho^2)}. \]
Vamos agora implementar o amostrador de Gibbs:
# parâmetros da simulação:
# número de amostras
m=1000
# coeficiente de correlação
rho=0.8
# valores médios:
mux=1
muy=2
# desvios padrão:
sigx=0.2
sigy=0.5
# inicialização de vetores:
x=rep(0,m)
y=rep(0,m)
#amostragem
# vamos começar de um valor arbitrário de y entre 0 e 1:
y[1]=runif(1)
for(i in 2:m){
# P(x|Y):
mu1=mux-rho*sigx*(y[i-1]-muy)/sigy
sig1=sigx/(1-rho^2)
x[i]=rnorm(1,mean=mu1,sd=sig1)
# P(y|x):
mu2=muy-rho*sigy*(x[i]-mux)/sigx
sig2=sigy/(1-rho^2)
y[i]=rnorm(1,mean=mu2,sd=sig2)
}
# vizualização:
plot(x,y)
points(mux,muy,col='red',pch=19,cex=1.5)
8 Ajuste de uma distribuição gaussiana
Vamos ver como podemos ajustar uma gaussiana a um conjunto de pontos:
# vamos simular 100 pontos com uma distribuição gaussiana de média 1.5 e desvio padrão 0.5:
x = rnorm(100, mean = 1.5, sd = 0.5)
# vamos carregar a biblioteca MASS:
library(MASS)
# vamos usar a função fitdistr do pacote MASS
gaussiana = fitdistr(x,'normal')
gaussiana## mean sd
## 1.40784244 0.56129840
## (0.05612984) (0.03968979)# parâmetros da gaussiana ajustada:
mu = gaussiana$estimate[[1]]
sig = gaussiana$estimate[[2]]
nn = gaussiana$n
# para ajustar a gaussiana:
xx = seq(-3,4,0.1)
h = hist(x, xlab='x', ylab='N')
h## $breaks
## [1] -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0
##
## $counts
## [1] 1 3 19 35 26 14 2
##
## $density
## [1] 0.02 0.06 0.38 0.70 0.52 0.28 0.04
##
## $mids
## [1] -0.25 0.25 0.75 1.25 1.75 2.25 2.75
##
## $xname
## [1] "x"
##
## $equidist
## [1] TRUE
##
## attr(,"class")
## [1] "histogram"lines(xx,dnorm(xx,mean=mu,sd=sig)*max(h$counts),lwd = 2)