In [1]:
# rm(list=ls())
options(OutDec = ",")
#==============================================================================
# Exemplo:
#
# Integracao Monte Carlo
#
#==============================================================================
set.seed(12345) # Semente
In [2]:
#==============================================================================
# Funcao
#==============================================================================
hx <- function(x){
(cos(50*x)+sin(20*x))^2
}
In [3]:
#==============================================================================
# Grafico da funcao a ser integrada
#==============================================================================
xseq <- seq(0,1,length=10001)
yseq <- hx(xseq)
par(mfrow=c(1,1),lwd=2,cex.lab=1.5,cex.axis=1.5,lab=c(10,10,0),
mar=c(4.5,5,2,1),bty="n")
plot(xseq,yseq,type="l",col="blue",lwd=3,xlim=c(0,1),ylim=c(0,4),
xlab=expression(x),ylab=expression(h(x)))
In [4]:
#==============================================================================
# Para n = 100 # Cem
#==============================================================================
n <- 100
u <- runif(n)
y <- hx(u)
print("===============================")
print(mean(y))
print("===============================")
print(sqrt(var(y)/n)) # Erro de Monte Carlo
[1] "===============================" [1] 0,8469294 [1] "===============================" [1] 0,09669878
In [5]:
#==============================================================================
# Para n = 1000 # Mil
#==============================================================================
n <- 1000
u <- runif(n)
y <- hx(u)
print("===============================")
print(mean(y))
print("===============================")
print(sqrt(var(y)/n)) # Erro de Monte Carlo
[1] "===============================" [1] 0,9254422 [1] "===============================" [1] 0,03230579
In [6]:
#==============================================================================
# Para n = 10000 # Dez mil
#==============================================================================
n <- 10000
u <- runif(n)
y <- hx(u)
print("===============================")
print(mean(y))
print("===============================")
print(sqrt(var(y)/n)) # Erro de Monte Carlo
[1] "===============================" [1] 0,9788495 [1] "===============================" [1] 0,01055842
In [7]:
#==============================================================================
# Para n = 100000 # Cem mil
#==============================================================================
n <- 100000
u <- runif(n)
y <- hx(u)
print("===============================")
print(mean(y))
print("===============================")
print(sqrt(var(y)/n)) # Erro de Monte Carlo
[1] "===============================" [1] 0,9672548 [1] "===============================" [1] 0,003306771
In [8]:
#==============================================================================
# Para n = 1000000 # Um milhao
#==============================================================================
n <- 1000000
u <- runif(n)
y <- hx(u)
print("===============================")
print(mean(y))
print("===============================")
print(sqrt(var(y)/n)) # Erro de Monte Carlo
[1] "===============================" [1] 0,964885 [1] "===============================" [1] 0,001045398
In [9]:
#==============================================================================
# Como exercicio, utilize a regra do trapezio, a de Simpson e a de Riemann
# para calcular a integral acima. Utilize diferentes numeros de subintervalos
# (e.g., 100, 1000, etc.)
#==============================================================================
In [10]:
# rm(list=ls())
#==============================================================================
# graphics.off()
#==============================================================================
# Fim
#==============================================================================