Estatística Computacional II - Lista 2
Aluno: Rafael Cabral Fernandez
Introdução
Este documento refere-se a lista 3 de exercícios proposta pelo professor Gustavo Rocha, da cadeira de Estatística Computacional 2, como cômputo de ensino da disciplina.
1
- Escreva um algoritmo para estimar \(P(X>2)\) onde \(X \sim Cauchy\), cuja densidade é \(p(x) = \cfrac{1}{\pi(1+x^2)}; x\in \mathbb{R}\). Considere integração de Monte Carlo e a densidade de amostragem por importância \(h(y) = \cfrac{2}{y^2}, y>2\)
Fx <- function(x){
2/(1-x)
}
Gx <- function(y){
2/( pi*(1+y^2)*(1/y^2) )
}
sample <- runif(1000000)
GXi <- Gx(sample)
mean(GXi)## [1] 0.13663152
Seja \(X\) um vetor aleatório de dimensão 3 com distribuição multinomial com parâmetros \(n = 5\), \(p_1= 0,5\), \(p_2= 0,2\) e \(p_3= 0,3\).
matriz <- matrix()
a<-vector()
b<-vector()
c<-vector()
for (i in 1:30000){
a<-rbinom(i,5,0.5)
b<-rbinom(i,5-a[i],2/5)
c<-rbinom(i,5-a[i]-b[i],3/5)
}
matriz <- rbind(a,b,c)
matrizfinal <- matriz[,which(colSums(matriz)==5)]
matrizfinal <- t((as.data.frame(matrizfinal)))
count(matrizfinal)## x.a x.b x.c freq
## 1 0 3 2 20
## 2 1 2 2 470
## 3 1 3 1 148
## 4 2 1 2 1475
## 5 2 2 1 1331
## 6 2 3 0 95
## 7 3 0 2 754
## 8 3 1 1 1889
## 9 3 2 0 450
## 10 4 0 1 472
## 11 4 1 0 328
## 12 5 0 0 353
Seja \(X\) uma variável aleatória contínua com função de densidade \(f(x)\) dada por \(f(x) = e^{-x}(12,5x^2e^{-4x} + 0,8); x>0\)
z1 <- rgamma(10000,3,5)
z2 <- rgamma(10000,1,1)
X = vector()
x1 = z1
x2 = z2
u = runif(10000,1)
if (u <= 0.2) {
X = z2
}else{
X = z1
}
hist(X,breaks = 50, probability = T)
4
Seja \(X \sim beta(1,4)\)
beta <- function(x){
1 - (1-x)^(1/4)
}
u<-runif(10000)
x<-beta(u)
eixox<-seq(0,1,l=10000)
eixoy<-dbeta(eixox,1,4)
hist(x,breaks = 50, probability = T)
lines(eixox,eixoy)
i<-1
j<-1
x<-c()
while (i<10000) {
U<-runif(1)
y<-runif(1)
r<-dbeta(y,1,4)/dunif(y,0,1)
c<-4
if(U<=r/c){
x[i]<-y
i=i+1}
j<-j+1
}
j## [1] 39864eixox<-seq(0.01,0.99,l=10000)
eixoy<-dbeta(eixox,1,4)
hist(x, breaks = 50, probability = T)
lines(eixox,eixoy)
#Taxa de Aceitação da amostra
10000/j## [1] 0.2508529#Taxa de Aceitação Teórica
1/c## [1] 0.25O valor de aceitação teórica é de 0.25 enquanto a taxa de aceitação está em 0.2498. Uma boa aproximação.
5
Estime a probabilidade \(P(0<Z<1)\), onde \(Z \sim N(0,1)\) para $n=10,100,10.00 e 10.000
- Método da aproximação retangular.
n <- 10
i <- seq(1:n)
h <- 1/n
ret10<-h*sum(dnorm(0+(2*i-1)/2*h))
n <- 100
i <- seq(1:n)
h <- 1/n
ret100<-h*sum(dnorm(0+(2*i-1)/2*h))
n <- 1000
i <- seq(1:n)
h <- 1/n
ret1000<-h*sum(dnorm(0+(2*i-1)/2*h))
n <- 10000
i <- seq(1:n)
h <- 1/n
ret10000<-h*sum(dnorm(0+(2*i-1)/2*h))
ret10## [1] 0.3414456ret100## [1] 0.3413458ret1000## [1] 0.3413448ret10000## [1] 0.3413447- Método trapezoidal.
a<-0
b<-1
n <- 10
i <- seq(1:n)
h <- (b-a)/n
trap10<-h*(dnorm(a)/2 + sum(dnorm(a+i*h)) + dnorm(b)/2)
n <- 100
i <- seq(1:n)
h <- (b-a)/n
trap100<-h*(dnorm(a)/2 + sum(dnorm(a+i*h)) + dnorm(b)/2)
n <- 1000
i <- seq(1:n)
h <- (b-a)/n
trap1000<-h*(dnorm(a)/2 + sum(dnorm(a+i*h)) + dnorm(b)/2)
n <- 10000
i <- seq(1:n)
h <- (b-a)/n
trap10000<-h*(dnorm(a)/2 + sum(dnorm(a+i*h)) + dnorm(b)/2)
trap10## [1] 0.3653401trap100## [1] 0.3437624trap1000## [1] 0.3415867trap10000## [1] 0.3413689- Regra de Simpson.
a<-0
b<-1
n <- 10
i1 <- seq(1:(n/2))
i2 <- seq(1:(n/2-1))
h <- (b-a)/n
simp10<-h/3*(dnorm(a) + 4*sum(dnorm(a+(2*i1-1)*h)) + 2*sum(dnorm(a+(2*i2-1)*h)) + dnorm(b))
n <- 100
i1 <- seq(1:(n/2))
i2 <- seq(1:(n/2-1))
h <- (b-a)/n
simp100<-h/3*(dnorm(a) + 4*sum(dnorm(a+(2*i1-1)*h)) + 2*sum(dnorm(a+(2*i2-1)*h)) + dnorm(b))
n <- 1000
i1 <- seq(1:(n/2))
i2 <- seq(1:(n/2-1))
h <- (b-a)/n
simp1000<-h/3*(dnorm(a) + 4*sum(dnorm(a+(2*i1-1)*h)) + 2*sum(dnorm(a+(2*i2-1)*h)) + dnorm(b))
n <- 10000
i1 <- seq(1:(n/2))
i2 <- seq(1:(n/2-1))
h <- (b-a)/n
simp10000<-h/3*(dnorm(a) + 4*sum(dnorm(a+(2*i1-1)*h)) + 2*sum(dnorm(a+(2*i2-1)*h)) + dnorm(b))
simp10## [1] 0.3453737simp100## [1] 0.3418559simp1000## [1] 0.3413969simp10000## [1] 0.34135- Método de Monte Carlo simples.
n<-10
mcs10<-mean(dnorm(runif(n)))
n<-100
mcs100<-mean(dnorm(runif(n)))
n<-1000
mcs1000<-mean(dnorm(runif(n)))
n<-10000
mcs10000<-mean(dnorm(runif(n)))
mcs10## [1] 0.3295192mcs100## [1] 0.3396191mcs1000## [1] 0.3425079mcs10000## [1] 0.341996- Método de Monte Carlo gerando valores de uma exponencial (1)
Indicadora<-function(x, min=0, max=1)
{
as.numeric(x >= min) * as.numeric(x <= max)
}
n<-10
y<-rexp(n,1)
mcexp10<-mean(Indicadora(y)*dnorm(y)/dexp(y))
n<-100
y<-rexp(n,1)
mcexp100<-mean(Indicadora(y)*dnorm(y)/dexp(y))
n<-1000
y<-rexp(n,1)
mcexp1000<-mean(Indicadora(y)*dnorm(y)/dexp(y))
n<-10000
y<-rexp(n,1)
mcexp10000<-mean(Indicadora(y)*dnorm(y)/dexp(y))
mcexp10## [1] 0.2385696mcexp100## [1] 0.3307061mcexp1000## [1] 0.3343831mcexp10000## [1] 0.3408038- Método de Monte Carlo gerando valores da própria normal padrão.
Indicadora<-function(x, min=0, max=1)
{
as.numeric(x >= min) * as.numeric(x <= max)
}
n<-10
y<-rnorm(n)
mcnorm10<-mean(Indicadora(y))
n<-100
y<-rnorm(n)
mcnorm100<-mean(Indicadora(y))
n<-1000
y<-rnorm(n)
mcnorm1000<-mean(Indicadora(y))
n<-10000
y<-rnorm(n)
mcnorm10000<-mean(Indicadora(y))
mcnorm10## [1] 0.3mcnorm100## [1] 0.36mcnorm1000## [1] 0.332mcnorm10000## [1] 0.3454- Para cada item acima faça um gráfico do valor da probabilidade em função de n. Ordene os melhores métodos para cada n.
plot(c(ret10,ret100,ret1000,ret10000),ylab="Valor",xlab="Método",xaxt="n",type="b",col="red")
axis(1,at=1:4,labels=c("N = 10","N = 100","N = 1000","N = 10000"))
main=c("n = 10")
title(main = "Método do Retângulo")
c=pnorm(1)-pnorm(0)
abline(h=c)
plot(c(trap10,trap100,trap1000,trap10000),ylab="Valor",xlab="Método",xaxt="n",type="b",col="red")
axis(1,at=1:4,labels=c("N = 10","N = 100","N = 1000","N = 10000"))
main=c("Método do Trapézio")
title(main = "Método do Trapézido")
c=pnorm(1)-pnorm(0)
abline(h=c)
plot(c(simp10,simp100,simp1000,simp10000),ylab="Valor",xlab="Método",xaxt="n",type="b",col="red")
axis(1,at=1:4,labels=c("N = 10","N = 100","N = 1000","N = 10000"))
main=c("Método de Simpsons")
title(main = "Método de Simpsons")
c=pnorm(1)-pnorm(0)
abline(h=c)
plot(c(mcs10,mcs100,mcs1000,mcs10000),ylab="Valor",xlab="Método",xaxt="n",type="b",col="red")
axis(1,at=1:4,labels=c("N = 10","N = 100","N = 1000","N = 10000"))
main=c("n = 10")
title(main = "Monte Carlo Simples")
c=pnorm(1)-pnorm(0)
abline(h=c)
plot(c(mcexp10,mcexp100,mcexp1000,mcexp10000),ylab="Valor",xlab="Método",xaxt="n",type="b",col="red")
axis(1,at=1:4,labels=c("N = 10","N = 100","N = 1000","N = 10000"))
main=c("n = 10")
title(main = "Método Carlo por Exponencial")
c=pnorm(1)-pnorm(0)
abline(h=c)
plot(c(mcnorm10,mcnorm100,mcnorm1000,mcnorm10000),ylab="Valor",xlab="Método",xaxt="n",type="b",col="red")
axis(1,at=1:4,labels=c("N = 10","N = 100","N = 1000","N = 10000"))
main=c("n = 10")
title(main = "Monte Carlo por Normal")
c=pnorm(1)-pnorm(0)
abline(h=c)
ordena <- function(x){
which.min(abs(c-x))
}
n10 <- c(ret10,trap10,simp10,mcs10,mcexp10,mcnorm10)
n100 <- c(ret100,trap100,simp100,mcs100,mcexp100,mcnorm100)
n1000 <- c(ret1000,trap1000,simp1000,mcs1000,mcexp1000,mcnorm1000)
n10000 <- c(ret10000,trap10000,simp10000,mcs10000,mcexp10000,mcnorm10000)
ordena(n10)## [1] 1ordena(n100)## [1] 1ordena(n1000)## [1] 1ordena(n10000)## [1] 1As melhores performances foram obtidos do Método de Retângulos para n = 10 e, para todos os demais, o Método de Simpsons alcançou otimalidade.