Exemplo já visto

Disponivel em Introdução à teoria assintótica (Gauss M. Cordeiro): https://www.ime.usp.br/~abe/lista/pdftCtIOIA62A.pdf. Exata, Aproximação Normal e expansão de Edgeworth para a função de densidade da soma estocástica padronizada de 5 variáveis aleatórias exponenciais independentes e identicamente distribuidas (iid).

library(EQL) #-> Extended-Quasi-Likelihood-Function
## Warning: package 'EQL' was built under R version 3.5.3
## Loading required package: ttutils
# Setando parâmetros iniciais
n = 5
mu = n
sigma2 = n
rho3 = 2
rho4 = 6

# Função exata
fp = function(y,n){
  (sqrt(n)*(n+y*sqrt(n))^(n-1)*exp(-n-y*sqrt(n)))/factorial(n-1)
}

# Aproximação de Edgeworth
my_edgeworth = function(y, n, rho3, rho4){
  dnorm(y)*(1+(rho3*hermite(y,3)/(6*sqrt(n)))+((rho4*hermite(y,4))/(24*n))
            +((rho3*rho3*hermite(y,6))/(72*n)))
}

# Grid de valores exatos
y = seq(-2, 3, 0.5); y
##  [1] -2.0 -1.5 -1.0 -0.5  0.0  0.5  1.0  1.5  2.0  2.5  3.0
Exato = round(fp(y,5),4)

# Aproximação pela normal
Normal = round(dnorm(y,0,1), digits = 4)

# Aproximação pela expansão de Edgeworth
Edgeworth = round(my_edgeworth(y, 5, 2, 6),4)

data.frame(y, Exato, Normal, Edgeworth)
##       y  Exato Normal Edgeworth
## 1  -2.0 0.0043 0.0540    0.0178
## 2  -1.5 0.1319 0.1295    0.1473
## 3  -1.0 0.3428 0.2420    0.3329
## 4  -0.5 0.4361 0.3521    0.4335
## 5   0.0 0.3924 0.3989    0.3923
## 6   0.5 0.2875 0.3521    0.2891
## 7   1.0 0.1840 0.2420    0.1886
## 8   1.5 0.1068 0.1295    0.1038
## 9   2.0 0.0577 0.0540    0.0500
## 10  2.5 0.0295 0.0175    0.0280
## 11  3.0 0.0144 0.0044    0.0182

Comparação gráfica:

Pareto Type I

Densidade: \[f(x; x_m, \alpha) = \frac{\alpha x_m^{\alpha}}{x^{\alpha + 1}} \]

Para mais detalhes ver: https://en.wikipedia.org/wiki/Pareto_distribution.

library(EnvStats)
## Warning: package 'EnvStats' was built under R version 3.5.3
## 
## Attaching package: 'EnvStats'
## The following objects are masked from 'package:stats':
## 
##     predict, predict.lm
## The following object is masked from 'package:base':
## 
##     print.default
n = 10
a = 5
k = 2
mu = (a*k)/(a-1)
sigma2 = (a*k^2)/((a-1)^2*(a-2))
rho3 <- log((2*(1+a)*sqrt((a-2)/a))/(a-3))
rho4 <- log((6*(a^3+a^2-6*a-2))/(a*(a-3)*(a-4)))

m = 100000
obs = rpareto(m*n, location = k, shape = a)

x = seq(2.0001, 6, 0.1)
plot(x,dpareto(x, shape = a, location = k), type = "l", col = "cyan")

# Distribuição empirica para uma pareto
x <- qemp(p = seq(0, 1, len = 1000), obs = obs) 
y <- demp(x, obs) 

# Distribuição empírica da soma estocástica de paretos
x = matrix(0, m, n)
for(i in 1:m){
  x[i,] = round(rpareto(n, location = k, shape = a),4)
}

sn = apply(x, 1, sum)

# Padronizando a soma estocástica
snp = (sn-n*mu)/(sqrt(n*sigma2))

x <- qemp(p = seq(0, 1, len = 100), obs = snp) 
y <- demp(x, snp) 

plot(x, y, type = "n", ylim = c(0,0.6), xlim = c(-2,5), 
     xlab = "Value of Random Variable", 
     ylab = "Relative Frequency", 
     main = "Soma Estocástica Padronizada n = 10") 
lines(x, y, lwd = 2, col = "cyan")
lines(x, dnorm(x), col = "#C72AC7", lwd = 2)
lines(x, my_edgeworth(x, n, rho3 = rho3, rho4 = rho4))
legend("topright", legend=c("Distribuição empírica","Normal","Expansão Edgeworth"),
       col = c("cyan", "#C72AC7", "black"), lty=1, cex=0.8)

Pareto Type II (Lomax Distribution)

Densidade: \[f(x; \alpha, \lambda) = \frac{\alpha}{\lambda}[1 + \frac{x}{\lambda}]^{-\alpha}\] Para mais detalhes ver: https://en.wikipedia.org/wiki/Lomax_distribution.

library(Renext)
## Warning: package 'Renext' was built under R version 3.5.3
## Loading required package: evd
## Warning: package 'evd' was built under R version 3.5.3
# Setando parâmetros iniciais
n = 10
l = 2
a = 5

mu = l/(a-1)
sigma2 = (a*l^2)/((a-1)^2*(a-2))
rho3 <- log((2*(1+a)*sqrt((a-2)/a))/(a-3))
rho4 <- log((6*(a^3+a^2-6*a-2))/(a*(a-3)*(a-4)))

x = seq(0.0001, 6, 0.1)
plot(x,dlomax(x, scale = l, shape = a), type = "l", col = "cyan")

m = 100000
obs = rlomax(m*n, scale = l, shape = a)

# Distribuição empirica para uma pareto
x <- qemp(p = seq(0, 1, len = 1000), obs = obs) 
y <- demp(x, obs) 

# Distribuição empírica da soma estocástica de Lomax
x = matrix(0, m, n)
for(i in 1:m){
  x[i,] = rlomax(n, scale = l, shape = a)
}

sn = apply(x, 1, sum)

#Trabalhando com a soma estocástica padronizada empirica
snp = (sn-n*mu)/(sqrt(n*sigma2))

x <- qemp(p = seq(0, 1, len = 100), obs = snp) 
y <- demp(x, snp) 

plot(x, y, xlim = c(-2,5), ylim = c(0,0.5), type = "n", 
     xlab = "Value of Random Variable", 
     ylab = "Relative Frequency",
     main = "Soma Estocástica Padronizada n = 10") 
lines(x, y, lwd = 2, col = "cyan")
lines(x, dnorm(x), col = "#C72AC7", lwd = 2)
lines(x, my_edgeworth(x, n, rho3 = rho3, rho4 = rho4))
legend("topright", legend=c("Distribuição empírica","Normal","Expansão Edgeworth"),
       col = c("cyan", "#C72AC7", "black"), lty=1, cex=0.8)

Triangular

Mais detalhes ver em: https://en.wikipedia.org/wiki/Triangular_distribution.

n = 5
a = -1
b = 1
c = 0
mu = (a+b+c)/3
sigma2 = (a^2+b^2+c^2-a*b-a*c-b*c)/18
rho3 = (sqrt(2)*(a+b-2*c)*(2*a-b-c)*(a-2*b+c))/((5*(a^2+b^2+c^2-a*b-a*c-b*c))^(3/2))
rho4 = -(3/5)

library(triangle)
## Warning: package 'triangle' was built under R version 3.5.3
m = 100000
obs = rtriangle(m*n,a,b,c)

x <- qemp(p = seq(0, 1, len = 100), obs = obs) 
y <- demp(x, obs) 

# Distribuição empírica da soma estocástica de triangulares
x = matrix(0, m, n)
for(i in 1:m){
  x[i,] = round(rtriangle(n, a, b, c),4)
}

sn = apply(x, 1, sum)

x <- qemp(p = seq(0, 1, len = 100), obs = sn) 
y <- demp(x, sn) 

#Trabalhando com a soma estocástica padronizada empirica
snp = (sn-n*mu)/(sqrt(n*sigma2))

x <- qemp(p = seq(0, 1, len = 100), obs = snp) 
y <- demp(x, snp) 

plot(x, y, type = "n", 
     xlab = "Value of Random Variable", 
     ylab = "Relative Frequency",
     main = "Soma Estocástica Padronizada") 
lines(x, y, lwd = 2, col = "cyan")
lines(x, dnorm(x), col = "#C72AC7", lwd = 1)
lines(x, my_edgeworth(x, n, rho3 = rho3, rho4 = rho4))
legend("topright", legend=c("Distribuição empírica","Normal","Expansão Edgeworth"),
       col = c("cyan", "#C72AC7", "black"), lty=1, cex=0.8)