Introdução

Sabemos que a distribuição Normal é seguramente a mais importante distribuição contínua que podemos trabalhar. Importante para diversas práticas, podemos destacar o Teorema Central do Limite, o qual é um resultado fundamental em aplicações práticas e teóricas, que garante a média de um conjunto de dados (que tenha outra distribuição) converge para uma distribuição normal conforme o número de dados aumenta.

Inicialmente vamos olhar como é a “forma”" de uma Normal, simulando uma amostra com 5.000.000 de observações.

n<-5000000
x<-rnorm(n)
m<-mean(x)
std<-sqrt(var(x))
hist(x, density=20, breaks=20, prob=TRUE, 
     xlab="Variável", ylab="Densidade",
     main="Curva Normal")
curve(dnorm(x, mean=m, sd=std), 
      col="darkblue", lwd=2, add=TRUE, yaxt="n")

O processo que buscamos consiste em estimar a integral da função densidade de probabilidade da Normal, que é dada por:

\[f(x)=\frac{1}{\sigma \sqrt{2\pi }}e^{-\frac{(x-\mu )^2}{2\sigma ^2}}\]

Devemos relembrar que a função de densidade de probabilidade descreve a probabilidade relativa de uma variável aleatória tomar um certo valor.

Essa probabilidade é dada pela integral da densidade dessa variável sobre tal faixa. A integral é um instrumento criado originalmente para determinar a área sob uma curva no plano cartesiano.

Desta forma, vamos mostrar como calcular a integral de uma distribuição Normal(\(\mu\),\(\sigma\)), através de uma função com o nome \(\mathbf{integra(a,b,\mu,\sigma,n)}\), que possibilite a obtenção da integral definida entre dois pontos \(a\) e \(b\) quaisquer, de uma distribuição Normal com esperança matemática \(\mu\) e desvio padrão \(\sigma\). Vamos usar um método geométrico por soma de retângulos, similar a imagem acima, com \(n\) representando o número de retângulos na aproximação.

O Modelo

Assim, para calcular tal probabilidade em uma distribuição normal, vamos usar os seguintes comandos:

Primeiramente precisamos imputar a função densidade de probabilidade da Normal

normal<-function(x,m,std){
1/(sqrt(2*pi)*std)*exp(-(x-m)^2/(2*std^2)) ##a altura dependerá da f(x) da normal
}

Perceba que podemos utilizar a \(f(x)\) da Normal como altura de cada retângulo que será utilizada para completar a integração. Para termos a área desse retângulo, precisamos ainda da sua base, que será descrita genericamente como um certo delta: \[\delta=\frac{b-a}{n}\] onde o \(n\) será o número de retângulos no intervalo \([a,b]\), de forma que quanto maior o valor de n melhor será a aproximação do valor calculado pela integral.

Por sua vez, para conseguir sua área será o somatório da multiplicação da base de cada retângulo por cada altura (valor na função de distribuição), ou seja \[Área=\delta\displaystyle\sum_{i=0}^{n-1} f(x_i)\]

Para integrar a Normal, vamos determinar esse intervalo genério embaixo da curva de \(a\) até \(b\):

integra<-function(a,b,m,std,n){
  dlt<-(b-a)/n ## base
  vr<-0:(n-1) ##vetor que está "contando o número de retângulos"
  x<-a+vr*dlt #criando os x necessários para utilizar na função
  y<-normal(x,m,std) ##altura dos retângulos
  sum(y*dlt) ## somatória de todas bases*altura = Área
}

Comparando o desempenho com a função pnorm do R

Conforme orientação, queremos agora trabalhar com uma distribuição Normal com média \(\mu\) = 100 e desvio-padrão \(\sigma\) = 30.

Vamos ver como será calculada a área abaixo da curva entre os pontos \(a = 35\) e \(b = 72\), com n = 50.000.

integra(35,72,100,30,50000)
## [1] 0.1601911

Agora vamos comprara com o resultado da função padrão do R para esse tipo de operação:

pnorm(72,100,30)-pnorm(35,100,30)
## [1] 0.1601938

Perceba que os resultados ficam realmente próximos. Vamos agora reduzir a nossa amostra para 500 e ver se perdemos muita qualidade no nosso modelo:

integra(35,72,100,30,500)
## [1] 0.1599226

O resultado muda bastante. Por isso é importante ter um \(n\) bem grande para esse tipo de aproximação. Vamos testar nosso comando para mais outros pontos.

Inicialmente \(a = -\infty\) (ou muito grande) e \(b = 10\), com \(n = 500.000\).

integra(10000000,10,100,30,500000)
## [1] -0.996702

E vamos comprar com a função padrão do R

pnorm(10,100,30)-pnorm(10000000,100,30)
## [1] -0.9986501

Ainda estamos relativamente pertos. Vamos testar agora para a outra cauda, com \(a = 200\) e \(b = \infty\) (ou muito grande), e \(n = 5.000.000\).

integra(200,10000000,100,30,5000000)
## [1] 0.0004823714

Vamos comprar com a função padrão do R

pnorm(10000000,100,30)-pnorm(200,100,30)
## [1] 0.0004290603

Em suma, o modelo sugerido tem ótimo desempenho para aproximações do cálculo de probabilidade, mas é necessário \(n\) bastante elevado. Desta forma, até por conta da simplicidade, é sugerido usar a função padrão do R.