Comandos R para análises estatísticas

Capítulo 9: Noções de Simulação

9.1 Introdução

Neste capítulo, revisaremos alguns comandos de simulação, já apresentados nos capítulos 6 e 7.
Caso o leitor tenha interesse pelo assunto, uma boa referência para métodos de simulação mais sofisticados implementados no R é o livro: Introducing Monte Carlo Methods with R, Robert, C.P., and Casella, G. , Springer, 2010.

Exemplo 9.4

Como vimos no capítulo 7, o comando para gerar números aleatórios de uma distribuição Uniforme é o comando runif():

u<-runif(10,0,1) # gera 10 observações aleatórias de uma v.a. uniforme[0,1]
u
##  [1] 0.75967 0.20125 0.25881 0.99215 0.80735 0.55333 0.64641 0.31182
##  [9] 0.62182 0.32977

Ainda, como visto no capítulo 7, se quisermos o quantil de uma distribuição normal, por exemplo, a partir da sua f.d.a, o comando utilizado é o qnorm. Em outras palavras, a função qnorm é a função inversa da função pnorm, que utilizamos para calcular a f.d.a da distribuição normal. Desta forma, se quisermos simular uma distribuição normal com média 0 e variância 1, usando o resultado do Teorema 9.1, faremos:

u <- runif(1000,0,1) # gera 1000 observações aleatórias de uma v.a. uniforme[0,1]
x <- qnorm(u,mean=0, sd = 1) # Calcula os quantis para o vetor u simulado da uniforme.

par(mfrow=c(1,2))
hist(u, freq=FALSE, main="Histograma da amostra da \n distribuição Uniforme simulada", col="lightblue4",border="gray")
hist(x, freq=FALSE, main="Histograma da variável X simulada a partir \n do resultado do Teorema 9.1", col="lightblue4",border="gray")

Figura 9.A: Exemplo de simulação da distribuição normal padrão a partir da simulação de uma distribuição Uniforme[0,1].

Exemplo 9.6

Desta vez, diferente do exemplo da Figura 9.A, em que temos a função qnorm, não possuímos a função inversa da f.d.a. da variável X:

\[ F(x)=\left\{\begin{array}{cl} 0 &,& x<0 \\ x^2 &,& 0\le x<1 \\ 1 &,& x\ge 1 \end{array}\right. \]

No entando esta função possui uma inversão muito fácil de calcular, i.e., \(F^{-1}(u) = \sqrt{u}\). Utilizaremos este resultado para calcular a amostra da variável x a partir da simulação de uma v.a. Uniforme[0,1]:

u <- runif(1000,0,1) # gera 1000 observações aleatórias de uma v.a. uniforme[0,1]
x <- sqrt(u) # Calcula a f.d.a. inversa para o vetor u gerado.

par(mfrow=c(1,2))
hist(u, freq=FALSE, main="Histograma da amostra da \n distribuição Uniforme simulada", col="lightblue4",border="gray")
hist(x, freq=FALSE, main="Histograma da variável x simulada a partir \n do resultado do Teorema 9.1", col="lightblue4",border="gray")

Figura 9.B: Histogramas das variáveis uniforme e normal geradas(n=1000).

Note que a distribuição de X simulada é muito próxima da Figura 9.4, apresentada a seguir:

par(mfrow=c(1,2))
x<-seq(0,1.75,0.01) # Produzindo o vetor da abscissa  
f_x<-c(2*x[x<=1],0*x[x>1]) # Calculando f(x)
plot(x,f_x, col='darkred', ylab="f(x)",xaxt="n",type="l") #Grafico de f(x)
lines(x=c(0,1),y=c(2,2), col="gray",lty=2)

F_x<-c((x[x<=1])^2,rep(1,length(x))[x>1]) # Calculando F(x)
plot(x,F_x, col='darkblue', ylab="F(x)",xaxt="n",type="l") #Grafico de F(x)
lines(c(0,1),c(1,1), col="gray",lty=2) # Adicionando linas explicativas tracejadas
lines(c(1,1),c(0,1), col="gray",lty=2)
lines(c(0,.71),c(0.51,.51), col="gray",lty=2)
lines(c(.71,.71),c(0,.51), col="gray",lty=2) 
text(0,.51,"u=0.5", pos=4) # inserindo texto no grafico
text(0.71,0,"x=0.71")

Figura 9.4: F.d.p e f.d.a. da v.a. X do Exemplo 9.6

Exemplo 9.8

Para simular uma variável X bernoulli com probabilidade p, vamos escrever uma função que verifica se a variável uniforme gerada é maior ou menor que p. Se for menor, então x=0, caso contrário, x=1:

rbernoulli<-function(n=1,p=0.5){
  u<-runif(n,0,1) # simula n uniformes
  x<-1*(u<=p) #categoriza as uniformes simuladas de acordo com o valor de p desejado
  return(x) 
}
rbernoulli(p=0.52)
## [1] 0

Se quisermos então gerar 10 valores desta v.a., faremos:

u<-rbernoulli(n=10,p=0.52)
u
##  [1] 1 1 1 0 1 0 0 1 0 0

Como sabemos, a proporção de 1’s em uma amostra aleatória de uma variável \(X\sim Bernoulli(p)\) é um estimador \(\hat{p}\) para \(p\).

aa<-rbernoulli(n=1000, p=0.52) 
sum(aa)/1000
## [1] 0.5

Exemplo 9.9

Para simular uma distribuição binomial a partir de experimentos de Bernoulli, basta repetirmos o processo de geração das v.a. bernoulli 20 vezes, calculando a soma de sucessos em cada uma:

xbinomial<-numeric()
for( i in 1:20){
  xbinomial[i]<-sum(rbernoulli(n=10,p=0.52))
}
xbinomial
##  [1] 7 7 4 5 6 6 4 6 5 5 8 6 5 5 8 5 5 4 7 4

Para simular uma distribuição binomial diretamente, basta utilizar o comando rbinom como visto no capítulo 6:

rbinom(n=20, size = 10, p=0.52)
##  [1] 7 7 5 6 8 4 5 6 7 2 9 4 8 6 7 6 4 5 5 5

Exemplo 9.10

Novamente, utilizemos o Teorema 9.1 para construir uma variável exponencial a partir de uma v.a. uniforme e a f.d.a. inversa da f.d. exponencial:

u<-runif(n=5,0,1)
qexp(u,rate=1/2)
## [1] 2.79958 3.94946 0.70592 1.01043 1.69785

Como já vimos, uma outra forma é simplesmente utilizar a função rexp:

rexp(n=5,rate=1/2)
## [1] 5.29244 0.63873 1.40804 1.57235 0.12447

Exemplo 9.11

Para gerar uma v.a. normal pela transformação integral, basta utilizamos a função qnorm, como fizemos na Figura 9.A:

qnorm(0.230,mean=10, sd = .4) # Calcula os quantis para a normal com média 10 e sd = 0.4
## [1] 9.7045

Como o R é baseado no S-plus, os comando do S-plus para simulação de variáveis aleatórias são os mesmos. Repetiremos assim, os mesmos comandos do S-plus nos próximos exemplos. Tais comandos também já foram apresentados neste material nos capítulos 6 e 7.

Exemplo 9.13

Simular 10 valores de uma v.a. \(X\sim b(10,0.5)\) e 20 valores de uma v.a. \(Y\sim Poisson(1.7)\):

x<-rbinom(20,10,0.5)
x
##  [1] 7 4 8 9 8 7 6 5 5 8 5 4 5 5 5 8 5 4 2 7
y<-rpois(20,1.7)
y
##  [1] 2 1 2 2 0 0 0 0 1 1 0 3 4 3 2 0 0 1 2 1

Exemplo 9.14

x<-rbinom(20,10,0.5)
y<-rpois(20,1.7)
z<-runif(100,0,1)
b<-rbinom(15,1,0.7)

Os histogramas das v.a. simuladas são apresentados na Figura 9.8:

par(mfrow=c(2,2))
hist(x, col="lightblue4",border="white")
hist(y, col="lightblue4",border="white")
hist(z, col="lightblue4",border="white")
hist(b, col="lightblue4",border="white")

Figura 9.8: Histogramas das distribuições simuladas no Exemplo 9.14. R

Exemplo 9.16

z<-rnorm(500,0,1)
y<-rnorm(200,10,3)
t<-rt(500,35)
Exp<-rexp(500,2)
w<-rchisq(300,5)
f<-rf(500,10,12)

par(mfrow=c(3,2))
hist(z, col="lightblue4",border="white")
hist(y, col="lightblue4",border="white")
hist(t, col="lightblue4",border="white")
hist(Exp, col="lightblue4",border="white")
hist(w, col="lightblue4",border="white")
hist(f, col="lightblue4",border="white")

Figura 9.16: Histogramas de algumas distribuições geradas no Exemplo 9.16.

Exemplo 9.18

Vamos simular duas v.a. uniformes independentes \(U_1 \sim \mathcal{U}[0,1]\) e \(U_2 \sim \mathcal{U}[0,1]\) de tamanho n=200. e plotá-las no quadrado unitário, gerando assim a Figura 9.1:

u1<-runif(200,0,1)
u2<-runif(200,0,1)
plot(u1,u2, pch=20, col="darkblue")

Figura 9.1: Área de Uma Figura por Simulação.

Exemplo 9.19

Vimos no Capítulo 8 como simular distribuições Normais multivariadas (não necessariamente independentes) com o pacote mvtnorm. Para gerar normais multivariadas com componentes independentes a partir do método de Box-Müller, apresentamos o código a seguir:

u1<-runif(10000,0,1)
u2<-runif(10000,0,1)
P<-data.frame(u1,u2)
x<-qnorm(u1)
y<-qnorm(u2)

Note que x e y são variáveis com distribuição normal padrão. Se plotarmos o diagrama de dispersão dessas duas variáveis, e suas marginais, teremos:

scatter.marginal(x,y,breaks=30) # Função criada no Capitulo 8 para construir este tipo de gráfico

Figura 9.C: Diagrama de dispersão de X e Y e suas marginais, simulados com o método de Box-Müller.

Para plotar a Figura 9.11, utilizaremos a função scatterplot3d:

library(mvtnorm)
library(scatterplot3d)
scatterplot3d(x,y, dmvnorm(x=cbind(x,y)), highlight.3d=TRUE, pch=16, zlab="Z")

Figura 9.11: Distribuição normal padrão bidimensional gerada: superfície.

Para plotar as curvas de nível utilizaremos a função contour, como no capítulo anterior. Como esta função exige que os dados das coordenadas X e Y sejam ordenados e equidistantes, utilizaremos a função interp da biblioteca akima para organizar os dados da maneira adequada para a função contour:

library(akima) 
contour(interp(x,y, dmvnorm(x=cbind(x,y))))

Figura 9.11: Distribuição normal padrão bidimensional gerada: curvas de nível.

Capítulo Anterior | Introdução | Próximo Capítulo