In [1]:
# rm(list=ls())
options(OutDec = ",") 
#==============================================================================
# Exemplo:
#
# Integracao Monte Carlo por funcao de importancia
#
# Calculando P(Z>=4,5) para Z ~ N(0,4)
#
#==============================================================================
# Probabilidade exata: P(Z>=4,5)=P(Z<=-4,5)
#==============================================================================
print(pnorm(-4.5))

[1] 3,397673e-06
In [2]:
#==============================================================================
# Monte Carlo simples para varios tamanhos de amostra
#==============================================================================
set.seed(12345)  # Semente
In [3]:
#==============================================================================
# n = 1000       # MIL 10^3
#==============================================================================
n   <- 1000      # MIL 10^3
z   <- rnorm(n)
print(mean(z>=4.5))

[1] 0
In [4]:
#==============================================================================
# n = 10000      # DEZ MIL 10^4
#==============================================================================
n   <- 10000     # DEZ MIL 10^4
z   <- rnorm(n)
print(mean(z>=4.5))

[1] 0
In [5]:
#==============================================================================
# n = 100000     # CEM MIL 10^5
#==============================================================================
n   <- 100000    # CEM MIL 10^5  
z   <- rnorm(n)
print(mean(z>=4.5))

[1] 2e-05
In [6]:
#==============================================================================
# n = 1000000    # UM MILHAO 10^6
#==============================================================================
n   <- 1000000   # UM MILHAO 10^6
z   <- rnorm(n)
print(mean(z>=4.5))

[1] 4e-06
In [7]:
#==============================================================================
# n = 10000000   # DEZ MILHOES 10^7
#==============================================================================
n   <- 10000000  # DEZ MILHOES 10^7
z   <- rnorm(n)
print(mean(z>=4.5))

[1] 4,3e-06
In [8]:
#==============================================================================
# Monte Carlo via funcao de importancia
# Funcao de importancia t-Student com 1 grau de liberdade
#==============================================================================
gl     <- 1           # Grau de liberdade
n      <- 1000        # Tamanho da amostra
z      <- rt(n,gl)    # Distribuicao de importancia
zcauda <- z[z>4.5]    # Seleciona os valores > 4,5
print(length(zcauda)) # Contagem de valores > 4,5
prob   <- sum(dnorm(zcauda)/dt(zcauda,gl))/n
print(prob)

[1] 74
[1] 2,062304e-06
In [9]:
#==============================================================================
# Probabilidade exata: P(Z>=4,5)=P(Z<=-4,5)
#==============================================================================
print(pnorm(-4.5))

[1] 3,397673e-06
In [10]:
#==============================================================================
# Monte Carlo via funcao de importancia
# Funcao de importancia Laplace(1), (escala tau=1)
#==============================================================================
silence <- suppressPackageStartupMessages # Para omitir mensagens de alertas
silence(library("smoothmest"))

n      <- 1000        # Tamanho da amostra
z      <- rdoublex(n) # Distribuicao de importancia
zcauda <- z[z>4.5]    # Seleciona os valores > 4,5
print(length(zcauda)) # Contagem de valores > 4,5
prob   <- sum(dnorm(zcauda)/ddoublex(zcauda))/n
print(prob)

[1] 3
[1] 3,771955e-06
In [11]:
#==============================================================================
# Probabilidade exata: P(Z>=4,5)=P(Z<=-4,5)
#==============================================================================
print(pnorm(-4.5))

[1] 3,397673e-06
In [12]:
#==============================================================================
# Grafico para comparacao
#==============================================================================
xseq   <- seq(0,6,length=1001)
zseq   <- dnorm(xseq)
yseq   <- ddoublex(xseq)
wseq   <- dt(xseq,gl)

par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(8,5,5),
    mar=c(4.5,5,1,1),xpd=T,cex.main=2.0,bty="n")
plot(xseq,zseq,type="l",lwd=4,col="black",xlab=expression(x),
  ylab=expression(f(x)),ylim=c(0,0.6))
lines(xseq,yseq,lwd=4,col="red")
lines(xseq,wseq,lwd=4,col="blue")
lines(c(4.5,4.5),c(0,0.1),col="darkorange",lwd=3,lty=2)
legend("topleft",legend=c(
  "Normal padrão",
  "Laplace com escala igual a 1",
  "t-Student com 1 grau de liberdade"),
  col=c("black","red","blue"),lty=c(1,1,1),lwd=c(4,4,4),bty="n",cex=1.5)
No description has been provided for this image
In [13]:
#==============================================================================
# Grafico para comparacao
# Olhando mais para a cauda superior
#==============================================================================
xseq   <- seq(3,8,length=1001)
zseq   <- dnorm(xseq)
yseq   <- ddoublex(xseq)
wseq   <- dt(xseq,gl)

par(mfrow=c(1,1),lwd=2.0,cex.lab=1.5,cex.axis=1.5,lab=c(6,5,5),
    mar=c(4.5,5,1,1),xpd=T,cex.main=2.0,bty="n")
plot(xseq,zseq,type="l",lwd=4,col=1,xlab="x",
  ylab=expression(f(x)),ylim=c(0,0.04))
lines(xseq,yseq,lwd=4,col="red")
lines(xseq,wseq,lwd=4,col="blue")
lines(c(4.5,4.5),c(0,0.03),col="darkorange",lwd=3,lty=2)
legend("topleft",legend=c(
  "Normal padrão",
  "Laplace com escala igual a 1",
  "t-Student com 1 grau de liberdade"),
  col=c("black","red","blue"),lty=c(1,1,1),lwd=c(4,4,4),bty="n",cex=1.5)
No description has been provided for this image
In [14]:
# rm(list=ls()) 
#==============================================================================
# graphics.off()
#==============================================================================
# Fim
#==============================================================================