--- title: "distribuições de probabilidades" author: "lsj" date: "March 23, 2023" output: html_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` ### 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 > - exercícios --- --- ## **1 A DISTRIBUIÇÃO GAUSSIANA EM R** > distribuição gaussiana de média $\mu$ e variância $\sigma^2$: $$ P(x) \sim N(\mu,\sigma) = {1 \over \sqrt{2 \pi} \sigma} \exp\Big[-{(x-\mu)^2 \over 2 \sigma^2} \Big] $$ > funções úteis para a 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 função de distribuição de probabilidades (fdp ou pdf) > * 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: ```{r , echo=TRUE} # reprodutibilidade set.seed(1234) # densidade de probabilidades, distribuição cumulativa, percentil da distribuição cumulativa e amostra de números aleatórios 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 para uma distribuição gaussiana de média 0 e desvio padrão 1 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 100 números aleatórios para uma gaussiana de média 10 e desvio padrão 5 x = rnorm(100,mean=10,sd=5) hist(x,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 $\pm$0.5 magnitude da média? ```{r , echo=TRUE} # 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) # 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 # visualização: hist(x,main='simulação',col='skyblue2') # * 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: ```{r , echo=TRUE} 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): ```{r , echo=TRUE} dbinom(1, 15, 0.10) ``` > 2. Nessa mesma amostra, qual é a probabilidade de não mais que 2 objetos serem RL? ```{r , echo=TRUE} dbinom(0, 15, 0.10) + dbinom(1, 15, 0.10) + dbinom(2, 15, 0.10) # ou, usando a função cumulativa: pbinom(2, 15, 0.10) ``` --- --- ## **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: ```{r , echo=TRUE} # reprodutibilidade set.seed(321) # meia vida (em unidades arbitrárias) meiavida = 2 # taxa: decaimentos por unidade de tempo lambda = 1/meiavida lambda # número inicial de átomos n0 = 1000 # tempos de decaimento de cada átomo td = rexp(n0,rate=lambda) summary(td) 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 # valor teórico: exp(-lambda*10)=exp(-5)= exp(-5) # 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 ```{r , echo=TRUE} x = -1/2*log(1-0.9) x ``` --- --- ## **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} $$ ```{r , echo=TRUE} # 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? > a) a probabilidade de termos 19 quasares ou menos é ```{r , echo=TRUE} ppois(19, 15) #lower tail ``` > b) a probabilidade de 20 ou mais quasares por unidade de área é ```{r , echo=TRUE} ppois(20, 15, lower=FALSE) #upper tail ``` > que é equivalente a ```{r , echo=TRUE} 1-ppois(20, 15) ``` > > 2. 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: ```{r , echo=TRUE} dpois(4,lambda=3)*dpois(4,lambda=4) ``` > Vamos verificar fazendo uma simulação: ```{r , echo=TRUE} # 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] # simulação das contagens do contador 2 n2 = rpois(nsim,lambda=4) n2[1:10] # 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) # fração das simulações com n1 e n2 iguais a 4 sum(sel)/nsim # 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: ```{r , echo=TRUE} mean(n1+n2) ``` --- --- ## **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: ```{r , echo=TRUE} 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: ```{r , echo=TRUE} 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 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. ```{r , echo=TRUE} 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 distribuição de $z$ fazendo uma simulação: ```{r , echo=TRUE} par(mfrow=c(1,2)) m=1000 x=rnorm(m) y=rnorm(m) z=x/y hist(z,main="distribuição de z=x/y") zz = z[z > -5 & z < 5] hist(zz,main="distribuição de z=x/y",xlim=c(-5,5),breaks=100) # note os outliers! summary(z) ``` > Pode-se mostrar que $z$ obedece a uma distribuição de Cauchy: $h(z) = {1 \over \pi} {1 \over 1+z^2}$ ```{r , echo=TRUE} 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) # distribuição cumulativa y = pcauchy(x) plot(x,pnorm(x),main="distr. cumulativa",type="l",col='red') # gaussiana lines(x,y,col='blue') # cauchy # quantis x = seq(0,1,0.01) # 100 números uniformemente distribuídos entre 0 e 1 y = qcauchy(x) plot(x,y,main="vetor de quantis",type="l",col='blue') lines(x,qnorm(x),col='red') # amostra aleatória x=rcauchy(100) hist(x,main="distr. aleatória ",col='blue') ``` --- --- ## **7 Simulações do Teorema do Limite Central** ### 7.1 simulação de distribuições exponenciais > 1. Vamos simular n conjuntos de m pontos gerados com uma distribuição exponencial > 2. Para cada simulação calculamos a média dos pontos > 3. 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 ```{r , echo=TRUE} 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) sd(x) ``` > Agora um exemplo de programa para esta simulação: ```{r , echo=TRUE} 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: ```{r , echo=TRUE} m=1000000 x=rexp(m) y=exp(-x) # média e desvio padrão desses pontos: mean(x) sd(x) ``` > ### 7.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: ```{r , echo=TRUE} # 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) head(df) # figura ggplot(df, aes(n, sigma_m)) + geom_line() # para obter uma curva mais suave, repetimos o processo de amostragem muitas vezes calc_sigma_m_mean = function(n, x) { mean(replicate(1000, sd(sample(x, n, replace = TRUE))/sqrt(n))) } df = data.frame(n, sigma_m = sapply(n, calc_sigma_m_mean, x)) ggplot(df, aes(n, sigma_m)) + geom_line() ``` > ### 7.3 *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 é: ```{r , echo=TRUE} 1-pnorm(100,mean=0,sd=81.65) ``` > * 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 ```{r , echo=TRUE} # 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 # comparem com o valor teórico, ~0.22 ``` --- --- ## **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(): ```{r} #library("mvtnorm") library('mixtools') # 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) # 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)) # 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: ```{r , echo=TRUE} # 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 # simulação gbv = rmvnorm(N, media, Sigma) # alguns pontos simulados: head(gbv) # 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)) c(mu2,mean(Y)) c(s1,sd(X)) c(s2,sd(Y)) # coeficiente de correlação entre as duas variáveis: cor(X,Y) ``` > > 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)$: ```{r} # 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. > * 2. 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. > * 3. Considere a distribuição $P(x,y) = a(x-y)^2$, com $0 * * a. determine $a$ para que $P(x,y)$ seja uma distribuição de probabilidades (isto é, seja normalizada) > * * b. determine $P(x)$ e $P(y)$ > * 4. 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)$. ```{r , echo=TRUE} ``` ```{r , echo=TRUE} ``` --- ---