Introdução e contextualização do problema

Essa nota técnica tem o objetivo de apresentar e testar diferentes técnicas de integração por métodos númericos. A dificuldade de implementação sistêmica de métodos de integração formais (teóricos) acaba por tornar muito atrativos a obtenção de \(Y\,|\,Y=\int_{a}^{b}f(x)\).

Para tal apresentação será utilizada a função de distribuição normal com médio \(\mu\) e desvio-padrão \(\sigma\). Como técnicas de integração serão utilizadas a soma de retângulos (pela esquerda e direita), soma de trapézios (ou média da soma de retângulos pela esquerda e direita), aproximação quadrática (método de Simpson) e por simulação de monte carlo.

Objetivos da nota

Essa nota tem como objetivos 4 pontos principais:

  1. Apresentar uma função que utilize todas as técnicas supracitadas para cálculo da integral definida de uma função de distribuição normal;

  2. Verificar a eficiência da função para o método das soma de retângulos comparando com a função “pnorm” do R para uma \(N(\mu=100,\sigma=30)\) os casos abaixo:
    • \((a=35,b=72)\)
    • \((a=-\infty,b=10)\)
    • \((a=200,b=\infty)\)

  3. Qual o valor de \(n\) (número de subdivisões) para que hajam ao menos \(m\) casas decimais corretas de resultado.

Função para calcular \(\int_{a}^{b} f(x) dx\) quando \(f(x)=N~(\mu,\sigma)\)

A função apresentada abaixo tem a finalidade de, por meio de 5 metodologias distintas de integração numérica, obter o valor de:

\[\begin{aligned} \int_{a}^{b}\frac{1}{\sigma\sqrt(2\pi)}e^-\frac{1}{2}(\frac{x-\mu}{\sigma})^2dx \end{aligned}\]

Abaixo encontra-se a função e uma apresentação de seus argumentos:

#    lim_inf = limite inferior de integração (a)
#    lim_sup = limite superior de integração (b)
#          n = número de repartições da curva para aproximação
#esperanca_x = média da distribuição de x
#       dp_x = desvio padrão da distribuição de x
#     method = método de integração numérica utilizado, tal que:
#             RLE  = SOMA DOS RETANGULOS PELA ESQUERDA
#             RLD  = SOMA DOS RETANGULOS PELA DIREITA
#             TRAP = SOMA DOS TRAPÉZIOS OU MÉDIA DA SOMA DE RETÂNGULOS
#             QUAD = APROXIMAÇÃO QUADRÁTICA
#             MC   = MONTE CARLO

Integra <- function(lim_inf=-Inf, lim_sup=Inf, n, esperanca_x, dp_x, method="MC"){

est_R <- pnorm(lim_sup,mean=esperanca_x,sd=dp_x)-pnorm(lim_inf,mean=esperanca_x,sd=dp_x)
lim_inf_orig<-lim_inf
lim_sup_orig<-lim_sup
  
if(is.finite(lim_inf)==FALSE & is.finite(lim_sup)==FALSE){
cat("Sua área será de 1 pois está contemplando toda a distribuição.") 
}

#RESOLVENDO PROBLEMA DE CALCULO PARA LIMITES DE INTEGRAÇÃO NO INFINITO
else if(is.finite(lim_inf)==FALSE & lim_sup <= esperanca_x){
  lim_inf = lim_sup
  lim_sup = esperanca_x
}

else if (is.finite(lim_inf)==FALSE & lim_sup > esperanca_x){
  lim_inf=esperanca_x
}

else if (is.finite(lim_sup)==FALSE & lim_inf <= esperanca_x){
  lim_sup=esperanca_x
}

else if (is.finite(lim_sup)==FALSE & lim_inf > esperanca_x){
  lim_sup=lim_inf
  lim_inf=esperanca_x
}

l<-(lim_sup-lim_inf)/n
i<-0:n
x<-lim_inf+i*l
y<- 1/(dp_x*sqrt(2*pi))*exp(-(((x-esperanca_x)/dp_x)^2)/2)

#SOMA DOS RETANGULOS PELA ESQUERDA
if (method == "RLE") {
  area_RLE <- ((lim_sup-lim_inf)/n)*sum(y[1:length(y)-1])
  
  if(is.finite(lim_inf_orig)==FALSE & lim_sup_orig <= esperanca_x){ 
  area_RLE<- 0.5-area_RLE}
  else if(is.finite(lim_inf_orig)==FALSE & lim_sup_orig > esperanca_x){ 
  area_RLE<- 0.5+area_RLE}
  else if(is.finite(lim_sup_orig)==FALSE & lim_inf_orig <= esperanca_x){ 
  area_RLE<- 0.5+area_RLE}
  else if(is.finite(lim_sup_orig)==FALSE & lim_inf_orig > esperanca_x){ 
  area_RLE<- 0.5-area_RLE}
    
erro_RLE<-abs(area_RLE-est_R)/est_R*100

cat("Method=", method, "area=", area_RLE, "est R=", est_R,"erro %=", erro_RLE,"\n") 
invisible(area_RLE)
}

#SOMA DOS RETANGULOS PELA DIREITA
 else if (method == "RLD") {
  area_RLD <- ((lim_sup-lim_inf)/n)*sum(y[2:length(y)])
  
  if(is.finite(lim_inf_orig)==FALSE & lim_sup_orig <= esperanca_x){ 
  area_RLD<- 0.5-area_RLD}
  else if(is.finite(lim_inf_orig)==FALSE & lim_sup_orig > esperanca_x){ 
  area_RLD<- 0.5+area_RLD}
  else if(is.finite(lim_sup_orig)==FALSE & lim_inf_orig <= esperanca_x){ 
  area_RLD<- 0.5+area_RLD}
  else if(is.finite(lim_sup_orig)==FALSE & lim_inf_orig > esperanca_x){ 
  area_RLD<- 0.5-area_RLD}
  
  erro_RLD<-abs(area_RLD-est_R)/est_R*100
  
cat("Method=", method, "area=", area_RLD, "est R=", est_R,"erro %=", erro_RLD,"\n") 
invisible(area_RLD)
}

#SOMA DOS TRAPÉZIOS
else if (method == "TRAP") {
  aux1 <- c(1,rep(2,n-1),1)
  area_trap <- ((lim_sup-lim_inf)/(2*n))*sum(y*aux1)
  
  if(is.finite(lim_inf_orig)==FALSE & lim_sup_orig <= esperanca_x){ 
  area_trap<- 0.5-area_trap}
  else if(is.finite(lim_inf_orig)==FALSE & lim_sup_orig > esperanca_x){ 
  area_trap<- 0.5+area_trap}
  else if(is.finite(lim_sup_orig)==FALSE & lim_inf_orig <= esperanca_x){ 
  area_trap<- 0.5+area_trap}
  else if(is.finite(lim_sup_orig)==FALSE & lim_inf_orig > esperanca_x){ 
  area_trap<- 0.5-area_trap}
  
  erro_trap<-abs(area_trap-est_R)/est_R*100
  
cat("Method=", method, "area=", area_trap, "est R=", est_R,"erro %=", erro_trap,"\n") 
invisible(area_trap)
}

#APROXIMAÇÃO QUADRÁTICA
else if (method == "QUAD") {
  aux2 <- c(c(1,4),rep(c(2,4),(n-2)/2),4)
  area_quad <- ((lim_sup-lim_inf)/(3*n))*sum(y*aux2)
  
  if(is.finite(lim_inf_orig)==FALSE & lim_sup_orig <= esperanca_x){ 
  area_quad<- 0.5-area_quad}
  else if(is.finite(lim_inf_orig)==FALSE & lim_sup_orig > esperanca_x){ 
  area_quad<- 0.5+area_quad}
  else if(is.finite(lim_sup_orig)==FALSE & lim_inf_orig <= esperanca_x){ 
  area_quad<- 0.5+area_quad}
  else if(is.finite(lim_sup_orig)==FALSE & lim_inf_orig > esperanca_x){ 
  area_quad<- 0.5-area_quad}
  
  erro_quad<-abs(area_quad-est_R)/est_R*100
  
  cat("Method=", method, "area=", area_quad, "est R=", est_R,"erro %=", erro_quad,"\n")
  invisible(area_quad)
  }

#SIMULAÇÃO DE MONTE CARLO
else if (method == "MC") {
  area_mc <- ((lim_sup-lim_inf)/n)*sum(y)
  
  if(is.finite(lim_inf_orig)==FALSE & lim_sup_orig <= esperanca_x){ 
  area_mc<- 0.5-area_mc}
  else if(is.finite(lim_inf_orig)==FALSE & lim_sup_orig > esperanca_x){ 
  area_mc<- 0.5+area_mc}
  else if(is.finite(lim_sup_orig)==FALSE & lim_inf_orig <= esperanca_x){ 
  area_mc<- 0.5+area_mc}
  else if(is.finite(lim_sup_orig)==FALSE & lim_inf_orig > esperanca_x){ 
  area_mc<- 0.5-area_mc}
  
  erro_mc <- abs(area_mc-est_R)/est_R*100
  
  cat("Method=", method, "area=", area_mc, "est R=", est_R,"erro %=", erro_mc,"\n") 
  invisible(area_mc)
}

}

Eficiência do método

  1. Para uma \(f(x)\approx N(\mu=100,\sigma=30)\) e limites de integração \((a=35,b=72)\) o resultado, para um cojunto de \(n=(10,100,1000)\) é:
Integra (35,72,10, 100, 30, method="TRAP")
## Method= TRAP area= 0.1603946 est R= 0.1601938 erro %= 0.1253196
Integra (35,72,100, 100, 30, method="TRAP")
## Method= TRAP area= 0.1601958 est R= 0.1601938 erro %= 0.001251897
Integra (35,72,1000, 100, 30, method="TRAP")
## Method= TRAP area= 0.1601938 est R= 0.1601938 erro %= 1.251884e-05
  1. No caso das integrais com limites em infinito ocorre um problema de cálculo relacionado a não convergência numérica do valor da base dos retângulos. Para lidar com esse problemas é possível assumir a premissa de que \(\int_{-\infty}^{\mu}f(x)dx=\int_{\mu}^{\infty}f(x)dx=0.5\). Desse modo, pode-se calcular, simetricamente, o valor de \(\int_{-\infty}^{b}f(x)dx\) ou \(\int_{b}^{\infty}f(x)dx\) da seguinte forma:

Se \(b<\mu\): \[\begin{aligned} \int_{-\infty}^{b}f(x)dx=\int_{-\infty}^{\mu}f(x)dx-\int_{b}^{\mu}f(x)dx=0.5-\int_{b}^{\mu}f(x)dx \newline \int_{b}^{\infty}f(x)dx=\int_{\mu}^{\infty}f(x)dx+\int_{b}^{\mu}f(x)dx=0.5+\int_{b}^{\mu}f(x)dx \end{aligned}\]

Se \(b>\mu\)::

\[\begin{aligned} \int_{-\infty}^{b}f(x)dx=\int_{-\infty}^{\mu}f(x)dx+\int_{\mu}^{b}f(x)dx=0.5+\int_{\mu}^{b}f(x)dx \newline \int_{b}^{\infty}f(x)dx=\int_{\mu}^{\infty}f(x)dx-\int_{\mu}^{b}f(x)dx=0.5-\int_{\mu}^{b}f(x)dx \end{aligned}\]

A função “integra” apresentada contempla essas situações e realiza os cálculos de forma adequada. Para uma \(f(x)\approx N(\mu=100,\sigma=30)\) e limites de integração \((a=-\infty,b=10)\) o resultado, para um cojunto de \(n=(10,100,1000)\) é:

Integra (-Inf,10,10, 100, 30, method="TRAP")
## Method= TRAP area= 0.001448719 est R= 0.001349898 erro %= 7.320639
Integra (-Inf,10,100, 100, 30, method="TRAP")
## Method= TRAP area= 0.001350895 est R= 0.001349898 erro %= 0.07386307
Integra (-Inf,10,1000, 100, 30, method="TRAP")
## Method= TRAP area= 0.001349908 est R= 0.001349898 erro %= 0.0007386965

De modo similar podemos obter para quando os limites de integração são \((a=200,b=\infty)\):

Integra (200,Inf,10, 100, 30, method="TRAP")
## Method= TRAP area= 0.000475953 est R= 0.0004290603 erro %= 10.92916
Integra (200,Inf,100, 100, 30, method="TRAP")
## Method= TRAP area= 0.0004295363 est R= 0.0004290603 erro %= 0.1109263
Integra (200,Inf,1000, 100, 30, method="TRAP")
## Method= TRAP area= 0.0004290651 est R= 0.0004290603 erro %= 0.001109427

Valor de \(n\) necessário para que convergência em \(m\) casas decimais

O procedimento abaixo permite determinar o valor de \(n\) para um dado \(m\) correspondente à precisão em casas decimais que se deseja.No exemplo, é utilizado como limites de integração \((a=200,b=\infty)\) para uma \(f(x)\approx N(\mu=100,\sigma=30)\) e \(m=4\).

  m=4
  i=1
integra_norm <- Integra(lim_inf = 200, lim_sup = Inf, n=i, esperanca_x = 100, dp_x = 100, method = "TRAP")
## Method= TRAP area= 0.1795435 est R= 0.1586553 erro %= 13.16581
est_r <- pnorm(Inf,mean=100,sd=30)-pnorm(200,mean=100,sd=30)
 
while(round(integra_norm,m)==round(est_r,m)){
i=i+1
integra_norm <- Integra(lim_inf = 200, lim_sup = Inf, n=i, esperanca_x = 100, dp_x = 100, method = "TRAP")
est_r <- pnorm(Inf,mean=100,sd=30)-pnorm(200,mean=100,sd=30)
  }
  
  print(i)
## [1] 1