Aplicação do Método de Expansão de Edgeworth

Universidade Federal de Pernambuco
jean.cardoso@ufpe.br

---

1 Introdução

Segundo Small (2010), a ideia de expandir uma função em uma série para estudar suas propriedades já existe há muito tempo. Newton desenvolveu algumas das fórmulas padrões que usamos hoje e Euler nos deu algumas ferramentas poderosas para somar séries. Assim, as expansões em série são certamente mais antigas do que o próprio assunto da estatística. Portanto, não é surpreendente encontrar expansões em série sendo usadas como ferramentas analíticas em muitas áreas da estatística. Para a grande maioria, o assunto é quase sinônimo de teoria assintótica. No entanto, as expansões em série surgem em muitos contextos, tanto na teoria das probabilidades quanto na estatística inferencial. Se definirmos assintótica no sentido amplo como o estudo de funções ou processos quando certas variáveis assumem valores limitantes, então todas as expansões em série são essencialmente investigações assintóticas.

O objeto central deste trabalho é o estudo de expansões assintóticas para distribuições de probabilidade. A rota mais óbvia é por meio de uma expansão explícita da função de distribuição acumulada ou da função densidade de probabilidade. Normalmente, embora não necessariamente, padronizamos a variável aleatória primeiro, antes de fazer a expansão, de modo que o termo principal da expansão corresponda, por exemplo, a uma distribuição normal de média zero e variância um. Em muitos casos encontrar a integral ou soma de convoluções explicitamente e derivar qualquer expansão assintótica de forma explícita é inviável ou até mesmo impossível. Devido a esta dificuldade vários métodos de aproximação foram desenvolvidos. Neste trabalho apresentaremos a expansão de Edgeworth na seção 3. Na seção 4, uma aplicação será feita com o objetivo de mostrar a eficácia deste método. Tomaremos como base a distribuição exponencial ponderada, que será apresentada a seguir.

2 Uma Nova Classe de Distribuições Exponenciais Ponderadas

Nos últimos anos, várias distribuições estatísticas foram propostas para modelar dados de tempo de vida que exibem funções de taxa de falha não constantes. Uma dessas distribuições é a distribuição exponencial ponderada (WE) de dois parâmetros, introduzida por Gupta e Kundu (2009) baseada na ideia de Azzalini (1985) que consiste na inserção de um novo parâmetro de forma a uma distribuição base. Primeiramente, vamos fornecer a definição da distribuição WE.

Dizemos que uma variável aleatória \(X\) segue uma distribuição WE, com parâmetro de forma e escala \(\alpha > 0\) e \(\lambda > 0\), respectivamente. A função densidade de probabilidade de \(X\) é dada por:

\[\begin{align*}\label{E:wexp} f(x\mid\alpha,\lambda) = \frac{\alpha +1}{\alpha}\lambda e^{-\lambda x}(1 - e^{-\alpha \lambda x}); \quad x>0 \end{align*}\]

e \(0\) caso contrário.

As correspondentes função de distribuição acumulada e função de risco são dadas pelas seguintes expressões, respectivamente:

\[\begin{align*} F(x\mid\alpha,\lambda) &= \frac{\alpha -1}{\alpha}e^{-\lambda x}(\alpha + 1 - e^{-\alpha\lambda x}), \\ h(x\mid\alpha, \lambda) &= (\alpha + 1)\lambda\frac{1 - e^{-\alpha\lambda x}}{\alpha + 1 - e^{-\alpha\lambda x}}. \end{align*}\]

Uma vez que a FDP de \(X\) é sempre log-côncava (ver Figura 1) e \(f(x\mid\alpha,\lambda)\) desaparece quando \(x\rightarrow0\) ou \(x\rightarrow\infty\), a função densidade de probabilidade é sempre unimodal com moda igual a \(\frac{\ln(1+\alpha)}{\lambda\alpha}\). Note que se \(\alpha\rightarrow0\), a WE converge para a distribuição gama com parâmetro de forma igual a \(2\) e o parâmetro de escala igual a \(\lambda\). Quando \(\alpha\rightarrow\infty\), ela converge para a distribuição exponencial com parâmetro de escala \(\lambda\). Caso \(\alpha = 1\), a WE coincide com a distribuição exponencial generalizada, com parâmetro de forma igual a \(2\).

Figura 1: Comportamento da função densidade de probabilidade da distribuição WE.

Observe que a função de risco \(h(x\mid\alpha,\lambda)\) (ver Figura 2) é uma função crescente em \(x\) com \(h(0\mid\alpha,\lambda) = 0\) e \(h(\infty\mid\alpha,\lambda) = \lambda\).

Figura 2: Comportamento da função de risco da distribuição WE.

Ao analisarmos a Figura 1, é possível perceber que a inserção do parâmetro de forma \((\alpha)\) trouxe mais flexibilidade para a função densidade de probabilidade. No painel direito, podemos observar que, dado \(\alpha\) fixo, cada valor de \(\lambda\) nos fornece um comportamento diferente da distribuição WE. Por outro lado, perceba que no painel esquerdo, fixado \(\lambda\), cada \(\alpha\) nos fornece uma pequena variação de comportamento. Isso é útil na análise de dados, pois nos dá a possibilidade de ajustar de forma mais precisa o comportamento da distribuição, fazendo com que a modelagem consiga detectar pequenas nuances. Talvez a grande importância da inclusão desse novo parâmetro possa ser observada justamente na Figura 2, pois ao analisarmos o painel direito, vemos que, fixado \(\alpha\), para qualquer valor de \(\lambda\) o comportamento da distribuição é único. Observando agora o painel esquerdo da Figura 2 é possível notar que, fixado \(\lambda\), assim como no caso da função densidade de probabilidade, temos pequenas variações para cada novo valor de \(\alpha\). Mais uma vez, vale ressaltar que a inserção desse parâmetro trás versatilidade à função de risco, o que é muito útil na análise de dados de sobrevivência.

Agora vamos considerar os vários momentos da distribuição WE. Se \(X\) segue uma distribuição WE, então a função geradora de momentos pode ser obtida da seguinte forma:

\[\begin{align*} M(t) = \mathrm{E}[e^{tx}] = \frac{\lambda^2(\alpha + 1)}{(\alpha\lambda + \lambda - t)(\lambda - t)}, \quad \text{para } t < \lambda. \end{align*}\]

A função geradora de momentos existe na forma compacta e, além disso, todos os seus momentos existem. Muitas vezes trabalhar com a função geradora de cumulantes \(K(t) = \ln M(t)\) é mais conveniente, pois torna mais fácil encontrar a esperança, variância, assimetria e curtose da distribuição. Calculando \(K(t)\) temos:

\[\begin{align*} K(t) = \ln M(t) = \ln\left(\alpha+1\right)+2\,\ln\left(\lambda\right)-\ln\left(\alpha\,\lambda+\lambda-t\right)-\ln\left(\lambda-t\right). \end{align*}\]

As quatro primeiras derivadas de \(K(t)\) aplicadas em \(t = 0\) são:

\[\begin{align*} k_{1} &= \mu = \mathrm{E}(X) = \frac{\alpha+2}{\lambda\,\left(\alpha+1\right)},\qquad k_{2} = \sigma^2 = \mathrm{var}(X) = \frac{{\alpha}^{2}+2\,\alpha+2}{{\lambda}^{2}\left(\alpha+1\right)^{2}}, \\ k_{3} &= \frac{2\,{\alpha}^{3}+6\,{\alpha}^{2}+6\,\alpha+4}{{\lambda}^{3}\left(\alpha+1\right)^{3}}, \quad \text{ } \quad k_{4} = \frac{6\,{\alpha}^{4}+24\,{\alpha}^{3}+36\,{\alpha}^{2}+24\,\alpha+ 12}{{\lambda}^{4}\left( \alpha+1\right)^{4}} \end{align*}\]

e correspondem a esperança, variância, assimetria e curtose, respectivamente.

Note que a média e a variância são ambas funções decrescentes de \(\alpha\) e ambas decrescem do valor \(2\) para o valor \(1\). Embora não seja possível calcular explicitamente a mediana, é possível observar que

\[\begin{align*} \text{moda} = \frac{\ln(1+\alpha)}{\lambda\alpha} < \text{mediana} < \frac{\alpha + 2}{\lambda(\alpha + 1)} = \text{média}. \end{align*}\]

pois a distribuição WE é assimétrica.

3 Método de Expansão de Edgeworth

Francis Ysidro Edgeworth (1854-1926) foi um grande filósofo e economista político anglo-inglês. Parte significativa de sua pesquisa foi na área da economia, onde ele buscou uma estrutura matemática para a teoria utilitarista. Durante a década de \(1880\) ele fez grandes contribuições para os métodos estatísticos, ficando conhecido pelos estatísticos pela expansão que leva seu nome.

A expansão de Edgeworth é um aprimoramento da aproximação do limite central. Ela inclui expressões que envolvem os cumulantes da distribuição de interesse nos termos de ordem superior. Essa expansão é usada para aproximações das somas padronizadas de variáveis aleatórias univariadas independentes e identicamente distribuídas (iid).

Na teoria assintótica, essas expansões são importantes quando a integral de convolução referente à soma de variáveis aleatórias não pode ser calculada explicitamente ou quando há interesse em obter aproximações para probabilidades do tipo \(P(S_{n}\ge s)\) ou \(P(\overline{Y} \ge y)\).

Suponha que \(Y_{1}, \ldots , Y_{n}\) são realizações iid de \(Y\) e sejam \(S_{n} = \sum_{i=1}^{n}Y_{i}\), a soma estocástica e \(S_{n}^{*} = (S_{n} - n\mu)/(\sigma\sqrt{n})\), a soma estocástica padronizada. As seguintes fórmulas são as expansões de Edgeworth para as funções densidade de probabilidade e de distribuição acumulada da soma estocástica padronizada, respectivamente

\[\begin{align} f_{S_{n}^{*}}(y) &= \phi(y) \bigg[1 + \frac{\rho_{3}}{6\sqrt{n}}H_{3}(y) + \frac{\rho_{4}}{24n}H_{4}(y) + \frac{\rho_{3}^{2}}{72n}H_{6}(y)\bigg] + O(n^{\frac{-3}{2}}), \\ F_{S_{n}^{*}}(y) &= \Phi(y) - \phi(y)\bigg[\frac{\rho_{3}}{6\sqrt{n}}H_{2}(y) + \frac{\rho_{4}}{24n}H_{3}(y) + \frac{\rho_{3}^{2}}{72n}H_{5}(y)\bigg] + O(n^{\frac{-3}{2}}), \end{align}\]

onde \(\rho_{3} = \frac{k_{3}}{(k_{2})^{3/2}}\) e \(\rho_{4} = \frac{k_{4}}{(k_{2})^{2}}\) são, respectivamente, o terceiro e o quarto cumulante padronizado, \(\Phi(y)\) e \(\phi(y)\) são, respectivamente, a função de distribuição acumulada e a função densidade de probabilidade da distribuição normal padrão e \(H_{i}\) são os polinômios probabilísticos de Hermite, com \(i = {2, 3, 4, 5, 6}\). Para mais detalhes sobre a formulação das expansões de Edgeworth consultar Small (2010), Barndorff-Nielsen e Cox (1990), Kolassa e McCullagh (1990) e McCullagh (1987).

4 Aplicação

Para o início da nossa análise, primeiramente carregamos o pacote EQL. Apesar desse pacote já possuir uma função chamada edgeworth, que calcula a expansão de Edgeworth da média padronizada, da média e da soma de variáveis aleatórias iid, neste trabalho optamos por construir manualmente as expansões de Edgeworth. Nesse pacote também está implementada uma função chamada hermite que calcula como default os polinômios probabilísticos de Hermite. Essa função é usada em nossas implementações.

# Carregando o pacote
library(EQL)

## Loading required package: ttutils

Como apresentado na Seção 3, para calcularmos as expansões de Edgeworth é necessário o cálculo dos cumulantes padronizados. Para isso vamos implementar os quatro primeiros cumulantes da distribuição WE.

## Primeiros Cumulantes
# 1º cumulante
cumulant1   <-  function(alpha, lambda)
{
    t7 <- (2 + alpha)/lambda/(alpha + 1)
    return(t7)
}

# 2º cumulante
cumulant2   <-  function(alpha, lambda)
{
    t1  <- alpha^2
    t4  <- lambda^2
    t8  <- (alpha + 1)^2
    t10 <- 1 / t8 / t4 * (t1 + 2*alpha + 2)
    return(t10)
}

# 3º cumulante
cumulant3   <-  function(alpha, lambda)
{
    t1  <- alpha^2
    t7  <- lambda^2
    t11 <- alpha + 1
    t12 <- t11^2
    t15 <- 1/t12/t11/t7/lambda*(2*t1*alpha + 6*alpha + 6*t1 + 4)
    return(t15)
}

# 4º cumulante
cumulant4   <-  function(alpha, lambda)
{
    t1  <- alpha^2
    t2  <- t1^2
    t9  <- lambda^2
    t10 <- t9^2
    t14 <- (alpha + 1)^2
    t15 <- t14^2
    t17 <- 1/t15/t10*(24*t1*alpha + 24*alpha + 36*t1 + 6*t2 + 12)
    return(t17)
}

A seguir implementamos as expansões de Edgeworth para a densidade e para a acumulada da soma estocástica padronizada.

## Aproximação de Edgeworth
# Densidade
density_edgeworth   <-  function(y)
{
    phi     <-  dnorm(y)
    t2      <-  (rho3*hermite(y,3))/(6*sqrt(n))
    t4      <-  (rho4*hermite(y,4))/(24*n)
    t6      <-  ((rho3^2)*hermite(y,6))/(72*n)
    t8      <-  phi*(1+t2)
    t10     <-  phi*(1+t2+t4+t6)
    list(ONm = t8, ON = t10)
}

# Acumulada
cumulative_edgeworth    <-  function(y)
{
    PHI     <-  pnorm(y)
    phi     <-  dnorm(y)
    t2      <-  (rho3*hermite(y,2))/(6*sqrt(n))
    t4      <-  (rho4*hermite(y,3))/(24*n)
    t6      <-  ((rho3^2)*hermite(y,5))/(72*n)
    t8      <-  PHI - phi*(t2)
    t10     <-  PHI - phi*(t2+t4+t6)
    list(ONm = t8, ON = t10)
}

Fixamos os parâmetros de alpha e lambda como \(2\). Definimos como variáveis globais \(\rho_{3}\) e \(\rho_{4}\), que são o terceiro e o quarto cumulantes padronizados da distribuição WE e, também de forma global, iniciaremos com \(n\) sendo \(2\).

# Setando valores
n       <<- 2
alpha   <-  2
lambda  <-  2
rho3    <<- cumulant3(alpha, lambda)/(cumulant2(alpha, lambda)^(3/2))
rho4    <<- cumulant4(alpha, lambda)/(cumulant2(alpha, lambda)^2)

Nossos próximos passos são setar um grid de valores, calcular a densidade e a acumulada da distribuição normal padrão e das expansões de Edgeworth, via simulação gerar \(10000\) valores pseudo-aleatórios da distribuição WE para cada uma das \(n\) variáveis aleatórias, somar essas \(n\) variáveis e, por fim, padronizá-las.

# Grid de valores exatos
y <- seq(-2, 3, 0.1);

## Aproximação pela distribuição normal padrão
dNormal <- dnorm(y)
pNormal <- pnorm(y)

## Aproximação pela expansão de Edgeworth
dEdgeworthONm   <- density_edgeworth(y)$ONm
dEdgeworthON    <- density_edgeworth(y)$ON
pEdgeworthONm   <- cumulative_edgeworth(y)$ONm
pEdgeworthON    <- cumulative_edgeworth(y)$ON

## Aproximação Empírica
# Geração de valores aleatórios pela soma de duas exponenciais
rwexp       <- function(n,alpha,lambda)
{
    V       <- rexp(n, lambda);
    W       <- rexp(n, (alpha + 1)*lambda);
    resul   <- V + W
    return(resul)
}

# Setando diretório
set.seed(1212)
obs <- 1e4
x   <-  matrix(0,obs,n)
for(i in 1:n){
    x[,i]   <-  rwexp(obs, alpha, lambda)
}

# Soma estocástica
sn = apply(x,1,sum)

# Soma estocástica padronizada
snp = (sn-(n*cumulant1(alpha, lambda)))/sqrt(n*cumulant2(alpha, lambda))

Finalmente vamos carregar o pacote latex2exp para auxiliar a escrita em LaTeX das legendas. Usando o comando plot criamos os gráficos.

## Gráficos

# Função auxiliar
auxfunc <- function(x, Digits = 4, Width = 4)
{
  formatC(x, digits = Digits, width = Width, format = "f")
}

# Carregando o pacote 
library(latex2exp)

# Plotando as densidades
par(mfrow = c(1,2), mar = c(3.2, 3.2, 0.4, 0.4), cex = 1.4)
plot(density(snp), type = 'l', lwd = 2, col = 1, lty = 1, xaxt = "n", yaxt = "n", 
    xlab = "", ylab = "", main = "", ylim = c(0, 0.5), xlim = c(-2, 3))
densidade <- list(dNormal, dEdgeworthONm, dEdgeworthON)
aux <- c(2,3,4)
for(i in 1:3){
    lines(y, densidade[[i]], lwd = 2, col = 1, lty = aux[i])
}
axis(2, seq(0, 0.5, length.out = 5), auxfunc(seq(0, 0.5, length.out = 5), 2))
axis(1, seq(-2, 3, length.out = 5), auxfunc(seq(-2, 3, length.out = 5), 2))
mtext("y", side = 1, line = 2.0, cex = 1.5)
mtext("Funções Densidade de Probabilidade", side = 2, line =2.0, cex = 1.5)
legend("topright", lty = 1:4, bty = "n", pch = c(20, 20, 20, 20), 
    legend = c("Empírica", "Normal", TeX("EW $O(n^{\\frac{-1}{2}}$)"), 
        TeX("EW $O(n^{-1}$)")), cex = 1, inset = 0.02, lwd = 2)

# Plotando as acumuladas
plot(ecdf(snp), lwd = 2, col = 1, lty = 1, xaxt = "n", yaxt = "n", 
    xlab = "", ylab = "", main = "", ylim = c(0,1), xlim = c(-2, 3))
acumulada <- list(pNormal, pEdgeworthONm, pEdgeworthON)
aux <- c(2,3,4)
for(i in 1:3){
    lines(y, acumulada[[i]], lwd = 2, col = 1, lty = aux[i])
}
axis(2, seq(0, 1, length.out = 5), auxfunc(seq(0, 1, length.out = 5), 2))
axis(1, seq(-2, 3, length.out = 5), auxfunc(seq(-2, 3, length.out = 5), 2))
mtext("y", side = 1, line = 2.0, cex = 1.5)
mtext("Funções de Distribuição Acumulada", side = 2, line =2.0, cex = 1.5)
legend("bottomright", lty = 1:4, bty = "n", pch = c(20, 20, 20, 20), 
    legend = c("Empírica", "Normal", TeX("EW $O(n^{\\frac{-1}{2}}$)"), 
        TeX("EW $O(n^{-1}$)")), cex = 1, inset = 0.03, lwd = 2)

Figura 3: Aproximações de Edgeworth (n = 2), distribuição normal padrão e curva empírica.

As Figuras 4, 5 e 6 são construídas usando o mesmo código mostrado anteriormente, alterando apenas \(n\) para os valores \(5\), \(10\) e \(20\).

Figura 4: Aproximações de Edgeworth (n = 5), distribuição normal padrão e curva empírica.

Figura 5: Aproximações de Edgeworth (n = 10), distribuição normal padrão e curva empírica.

Figura 6: Aproximações de Edgeworth (n = 20), distribuição normal padrão e curva empírica.

Analisando o painel esquerdo da Figura 3, é possível observar que a EW \(O(n^{-1})\) se aproxima mais da curva empírica do que a EW \(O(n^{-\frac{1}{2}})\). Entretanto, ambas estão relativamente mais próximas da curva empírica do que a distribuição normal. No painel direito observamos um comportamento semelhante. Note que no painel esquerdo da Figura 4, as curvas EW \(O(n^{-\frac{1}{2}})\), EW \(O(n^{-1})\) e empírica estão quase sobrepostas, mas já começam a ficar mais próximas da densidade da normal. O mesmo ocorre para o painel direito. Na Figura 5, podemos ver que as curvas estão se aproximando cada vez mais da distribuição normal. Podemos verificar na Figura 6 que, a medida que estamos aumentando o valor de \(n\), a curva normal padrão e a curva empírica estão muito próximas, não sendo necessário utilizar as expansões de Edgeworth.

5 Conclusão

Neste trabalho, apresentamos a nova classe de modelos exponenciais ponderados, criada acrescentando um novo parâmetro de forma a uma distribuição não simétrica. Mostramos os comportamentos das funções densidade de probabilidade e de risco. É possível observar as mudanças que o novo parâmetro traz à distribuição base, dando maior flexibilidade e tornando esse modelo interessante para aplicações na área de análise de sobrevivência. Algumas propriedades são apresentadas, entre elas a função geradora de momentos, a função geradora de cumulantes, a média, a variância, a assimetria e a curtose todas em sua forma analítica. Apresentamos de forma sucinta as expansões de Edgeworth para a soma estocástica padronizada de variáveis aleatórias iid. Comparamos a distribuição empírica e as aproximações de Edgewoth à distribuição normal. Foi possível observar que, já para o valor \(n=2\), a curva de EW \(O(n^{-1})\) se ajustou bem a distribuição empírica. Vale salientar que a partir de \(n=5\), as curvas de EW \(O(n^{-\frac{1}{2}})\) e EW \(O(n^{-1})\) coincidem e além disso se ajustam muito bem a curva empírica. Isso mostra que mesmo a expansão de Edgeworth de ordem \(O(n^{-\frac{1}{2}})\) já nos oferece um ajuste satisfatório para um baixo valor de \(n\). Podemos observar que a partir de \(n=20\), por exemplo, a distribuição normal padrão oferece um ajuste razoável para a curva empírica, por isso, não é mais necessário utilizar as expansões de Edgeworth.

6 Referências

Azzalini, A. (1985). A Class of Distributions Which Includes the Normal Ones. Scandinavian Journal of Statistics, 12(2), 171-178.

Barndorff-Nielsen, O. E. and Cox, D. R. (1990). Asymptotic Techniques for use in Statistics. Londres: Chapman and Hall.

Gupta, Rameshwar D. and Kundu, Debasis (2009). A new class of weighted exponential distributions. Statistics, 43:6, 621-634.

Kolassa, J. and McCullagh, P. (1990). Edgeworth expansions for discrete distributions. Ann. Statist., 18, 981–985.

McCullagh, P. (1987). Tensor methods in Statistics. Londres: Chapman and Hall.

R Core Team (2019). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.

Small, Christopher G. (2010). Expansions and Asymptotics for Statistics (Chapman & Hall CRC Monographs on Statistics & Applied Probability). Chapman and Hall/CRC.