Aplicação do Método de Expansão Ponto de Sela
Jean Carlos Cardoso
Universidade Federal de Pernambucojean.cardoso@ufpe.br
01/07/2021
1 Introdução
Com o avanço da estatística como ciência, os métodos modernos passaram a utilizar modelos estatísticos sofisticados, que requerem cálculos de probabilidades extremamente complexos. Desse modo, muitas vezes esses cálculos se tornam inviáveis ou intratáveis do ponto de vista analítico. As expansões ponto de sela são uma resposta adequada para essa nova situação.
Daniels (1954) foi pioneiro na aplicação da expansão ponto de sela na área da estatística. Mais tarde Lugannani e Rice (1980) desenvolveram uma fórmula alternativa que depende não apenas do ponto de sela mas também de uma singularidade na origem.
Neste trabalho temos como objetivo o estudo das expansões ponto de sela criadas por Daniels (1954, 1987) e Lugannani e Rice (1980). Este trabalho está disposto da seguinte forma. Na seção 2 é apresentada a distribuição que usaremos como base para nossas aplicações. Na seção 3 apresentaremos as expansões ponto de sela para a soma de variáveis estocásticas e para a média amostral. Além disso, exploramos duas fórmulas para a aproximação da distribuição acumulada, uma criada por Daniels (1987) e outra desenvolvida por Lugannani e Rice (1980). A seção 4 apresenta uma aplicação. O objetivo da mesma é demonstrar como as expansões ponto de sela são eficientes para aproximações. Por fim, a seção 5 traz as conclusões obtidas neste texto.
2 Distribuição Exponencial Ponderada
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 à \(\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 à \(2\) e o parâmetro de escala igual à \(\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 à \(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 traz 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 Ponto de Sela
Henry Ellis Daniels (1912-2000) foi presidente da Royal Statistical Society de 1974 a 1975 e uma figura chave no desenvolvimento da teoria, prática e ensino da estatística. Ele combinou excepcionais habilidades em matemática aplicada com uma compreensão pragmática da estatística. Seu artigo de 1954 introduziu o método do ponto de sela na estatística (SMALL, 2010).
As expansões ponto de sela são muito importantes na teoria assintótica, pois aproximam de forma precisa as funções densidades e de distribuições, sendo facilmente deduzida de sua função geradora de cumulantes. Também utilizaremos uma fórmula alternativa proposta por Lugannani e Rice (1980) para a aproximação da função de distribuição.
Suponha que \(Y_{1}, \ldots , Y_{n}\) são realizações iid de \(Y\) e seja \(S_{n} = \sum_{i=1}^{n}Y_{i}\) a soma estocástica das variáveis aleatórias contínuas.
A fórmula da expansão ponto de sela para a aproximação da função densidade da soma estocástica é dada por
\[\begin{align} f_{S_{n}}(s) = \frac{\exp\left\{n\,K(\widehat{t}) - s\,\widehat{t}\right\}}{\sqrt{2\,n\,\pi\,K^{(2)}(\widehat{t})}}\left\{1 + M(\widehat{t}) + O(n^{-2})\right\}, \end{align}\]
onde \(M(t)\) é um termo de ordem \(n^{-1}\) dado por
\[\begin{align} M(t) = \frac{3\,\rho_{4}(t) - 5\,\rho_{3}(t)^{2}}{24\,n}, \end{align}\]
sendo \(\rho_{j}(t) = K^{(j)}(t)/K^{(2)}(t)^{j/2}\) para \(j=3\), \(4\) e \(K^{(j)}(t)=d^{j}K(t)/d\,t^{j}\). Assim \(\rho_{3}(t)\) e \(\rho_{4}(t)\) são os cumulantes padronizados que medem a assimetria e a curtose da distribuição base. Os termos \(K(t)\) e \(K^{2}(t)\) são, respectivamente, a função geratriz de cumulantes e sua derivada segunda. Por fim, \(\widehat{t}\) deve satisfazer a equação \(K'(\widehat{t}) = s/n\).
A expansão para a função densidade da média amostral \(\overline{Y}_{n} = S_{n}/n\) é dada por
\[\begin{align} f_{\overline{Y}_{n}}(y) = \left\{\frac{n}{2\,\pi\,K^{(2)}(\widehat{t})}\right\}^{1/2} \exp\left\{n\,K(\widehat{t})-n\,\widehat{t}\,y \right\}\left\{1 + M(\widehat{t}) + O(n^{-2})\right\}. \end{align}\]
Usualmente, o interesse maior em inferência reside em obter aproximações precisas para probabilidades do tipo \(P(S_{n}\geq s)\) ou \(P(\overline{Y}_{n} \geq y)\) de uma amostra iid de \(n\) observações. Daniels (1987) exibiu a expansão ponto de sela para \(P(S_{n}\geq s)\) até termos de ordem \(O(n^{-1})\) nos seguintes casos:
- Caso 1: Quando \(s>n\,E(Y)\), ou seja, \(\widehat{t}>0\) temos que
\[\begin{align} P(S_{n}\geq s) &= \exp\left\{n\,K(\widehat{t})-s\,\widehat{t}+\frac{\widehat{v}^{2}}{2}\right\}\times \\ &\times\Bigg\{\left[1-\Phi(\widehat{v})\right]\left[1-\frac{\rho_{3}(\widehat{t})\,\widehat{v}^{3}}{6\,\sqrt{n}}+\frac{1}{n}\left(\frac{\rho_{4}(\widehat{t})\,\widehat{v}^{4}}{24}+\frac{\rho_{3}(\widehat{t})^{2}\,\widehat{v}^{6}}{72}\right)\right] + \\ &+ \phi(\widehat{v})\left[\frac{\rho_{3}(\widehat{t})(\widehat{v}^{2}-1)}{6\,\sqrt{n}}-\frac{1}{n}\left(\frac{\rho_{4}(\widehat{t})(\widehat{v}^{3}-\widehat{v})}{24}+\frac{\rho_{3}(\widehat{t})^{2}(\widehat{v}^{5}-\widehat{v}^{3}+3\widehat{v})}{72}\right)\right] \Bigg\}, \end{align}\]
onde \(\widehat{v} = \widehat{t}\left[n\,K^{(2)}(\widehat{t})\right]^{1/2}\).
- Caso 2: Quando \(s<n\,E(Y)\), ou seja, \(\widehat{t}<0\), podemos obter \(P(S_{n}\geq s)\) até \(O(n^{-1})\) como
\[\begin{align} P(S_{n}\geq s) &= H(-\widehat{v}) + \exp\left\{n\,K(\widehat{t})-s\,\widehat{t}+\frac{\widehat{v}^{2}}{2}\right\} \times \\ &\times\Bigg\{\left[H(\widehat{v})-\Phi(\widehat{v})\right]\left[1-\frac{\rho_{3}(\widehat{t})\,\widehat{v}^{3}}{6\,\sqrt{n}}+\frac{1}{n}\left(\frac{\rho_{4}(\widehat{t})\,\widehat{v}^{4}}{24}+\frac{\rho_{3}(\widehat{t})^{2}\,\widehat{v}^{6}}{72}\right)\right]+ \\ &+ \phi(\widehat{v})\left[\frac{\rho_{3}(\widehat{t})(\widehat{v}^{2}-1)}{6\,\sqrt{n}}-\frac{1}{n}\left(\frac{\rho_{4}(\widehat{t})(\widehat{v}^{3}-\widehat{v})}{24}+\frac{\rho_{3}(\widehat{t})^{2}(\widehat{v}^{5}-\widehat{v}^{3}+3\widehat{v})}{72}\right)\right]\Bigg\}, \end{align}\]
onde \(H(w)=0\), \(1/2\) e \(1\) quando \(w<0\), \(w=0\) e \(w>0\), respectivamente.
Uma forma alternativa simples de obter \(P(S_{n}\leq s)\) é deduzida por Lugannani e Rice (1980)
\[\begin{align} P(S_{n}\leq s) = \Phi(\widehat{r}) + \left(\frac{1}{\widehat{r}} - \frac{1}{\widehat{v}}\right)\phi(\widehat{r}), \end{align}\]
onde \(\widehat{r} = \text{sinal}(\widehat{t})\left\{2\,n\left[\widehat{t}K^{'}(\widehat{t})-K(\widehat{t})\right]\right\}^{1/2}\), cujo erro é na ordem de \(O(n^{-1})\) uniformemente em \(s\).
A aproximação anterior é boa em quase todo o intervalo de variação de \(s\), exceto próximo ao ponto \(s=E(S_{n})\) ou \(r=0\), onde deve ser substituída pelo seu limite, quando \(r\rightarrow 0\), dada por
\[\begin{align} P(S_{n}\leq s) = \frac{1}{2} + \frac{\rho_{3}(\widehat{t})}{6\sqrt{2\,\pi\,n}}. \end{align}\]
Mais detalhes sobre as expansões acima mencionadas podem ser encontrados em Cordeiro (1999).
4 Aplicação
Neste trabalho optamos por utilizar o software R. Além disso, não será necessária a utilização de nenhum pacote específico visto, que implementaremos nossas próprias expansões ponto de sela.
Iniciamos nossa análise implementando a função geratriz de cumulantes bem como os quatro primeiros cumulantes.
## Primeiros Cumulantes
# Função Geratriz de Cumulantes
K <- function(t)
{
t1 <- lambda^2
t2 <- log(t1)
t4 <- log(lambda-t)
t5 <- -alpha-1
t10 <- log(0.1e1/(lambda*t5+t)*t5)
t11 <- t2-t4+t10
return(t11)
}
# Primeiro Cumulante
K1 <- function(t)
{
t1 <- alpha * lambda
t10 <- 1 / (lambda - t) / (t1 + lambda - t) * (t1 + 2 * lambda - 2 * t)
return(t10)
}
# Segundo Cumulante
K2 <- function(t)
{
t1 <- alpha^2
t4 <- lambda^2
t10 <- t^2
t14 <- (-lambda+t)^2
t20 <- ((-alpha-1)*lambda+t)^2
t22 <- 1/t20/t14*(t4*(t1+2*alpha+2)-2*t*(alpha+2)*lambda+2*t10)
return(t22)
}
# Terceiro Cumulante
K3 <- function(t)
{
t5 <- alpha^2
t7 <- lambda^2
t12 <- t^2
t15 <- -lambda+t
t16 <- t15^2
t21 <- (-alpha-1)*lambda+t
t22 <- t21^2
t28 <- -4/t22/t21/t16/t15*(t7*(t5+alpha+1)-t*(alpha+2)*lambda+t12)*
((-alpha/2-1)*lambda+t)
return(t28)
}
# Quarto Cumulante
K4 <- function(t)
{
t1 <- alpha^2
t2 <- t1^2
t9 <- lambda^2
t10 <- t9^2
t12 <- alpha+2
t19 <- t^2
t29 <- t19^2
t33 <- (-lambda+t)^2
t34 <- t33^2
t40 <- ((-alpha-1)*lambda+t)^2
t41 <- t40^2
t43 <- 1/t41/t34*(t10*(24*t1*alpha+24*alpha+36*t1+6*t2+12)-
24*t9*lambda*(t1+alpha+1)*t12*t+36*t9*(t1+2*alpha+2)*t19-
24*lambda*t12*t19*t+12*t29)
return(t43)
}Como apresentado na Seção 3, precisamos criar nossas expansões ponto de sela. A seguir criamos a função densitysaddlepoint, meansaddlepoint, danielssaddlepoint e LRsaddlepoint, que são as correspondentes aproximações para a densidade, a média, a acumulada criada por Daniels (1987) e a acumulada criada por Lugannani e Rice (1980). Vale ressaltar que nossa distribuição tem \(\widehat{t}<0\). Logo, usaremos o caso 2 da distribuição acumulada de Daniels (1987). Todas as funções aqui implementadas tem erro na ordem de \(O(n^{-1})\).
## Aproximação Ponto de Sela
# Densidade da Soma Estocástica
densitysaddlepoint <- function(s)
{
t <- -(-alpha*lambda*s-2*lambda*s+sqrt(alpha^2*lambda^2*s^2+4*n^2)+2*n)/(2*s)
rho3 <- K3(t)/(K2(t)^(3/2))
rho4 <- K4(t)/(K2(t)^2)
Mhat <- (3*rho4-5*((rho3)^2))/(24*n)
final <- exp(n*K(t)-t*s)*sqrt(1/(2*pi*n*K2(t)))*(1+Mhat)
return(final)
}
# Densidade da Média
meansaddlepoint <- function(s)
{
t <- -(-alpha*lambda*s-2*lambda*s+sqrt(alpha^2*lambda^2*s^2+4*n^2)+2*n)/(2*s)
rho3 <- K3(t)/(K2(t)^(3/2))
rho4 <- K4(t)/(K2(t)^2)
Mhat <- (3*rho4-5*((rho3)^2))/(24*n)
aux0 <- sqrt(n/(2*pi*K2(t)))
aux1 <- exp(n*K(t)-t*s)
aux2 <- 1+Mhat
final <- aux0*aux1*aux2
return(final)
}
# Função Auxiliar para a Distribuição Acumulada (Daniels)
signH <- function(s)
{
i <- 0
aux <- NULL
for(i in 1:length(s)){
if(s[i]<0){
aux[i] <- 0
} else if(s[i]==0){
aux[i] <- 1/2
} else{
aux[i] <- 1
}
i = i+1
}
return(aux)
}
# Distribuição Acumulada_Daniels (1987) (\hat{t} < 0)
danielssaddlepoint <- function(s)
{
t <- -(-alpha*lambda*s-2*lambda*s+sqrt(alpha^2*lambda^2*s^2+4*n^2)+2*n)/(2*s)
rho3 <- K3(t)/(K2(t)^(3/2))
rho4 <- K4(t)/(K2(t)^2)
v <- t*sqrt(n*K2(t))
aux0 <- exp(n*K(t)-s*t+(v^2)/2)
aux1 <- signH(v) - pnorm(v)
aux2 <- 1 - (rho3*(v^3))/(6*sqrt(n)) + (rho4*(v^4))/(24*n) +
((rho3^2)*(v^6))/(72*n)
aux3 <- (rho3*(v^2 - 1))/(6*sqrt(n)) - (rho4*(v^3 - v))/(24*n) -
((rho3^2)*(v^5 - v^3 +3*v))/(72*n)
final <- signH(-v) + aux0*(aux1*aux2 + dnorm(v)*aux3)
return(1-final)
}
# Distribuição Acumulada_Lugannani e Rice (1980)
LRsaddlepoint <- function(s)
{
t <- -(-alpha*lambda*s-2*lambda*s+sqrt(alpha^2*lambda^2*s^2+4*n^2)+2*n)/(2*s)
q <- sign(t)*sqrt(2*(t*K1(t)-K(t)))
r <- sqrt(n)*q
v <- t*sqrt(n*K2(t))
final <- pnorm(r) + (1/(r) - 1/(v))*dnorm(r)
return(final)
}Fixamos, de forma global, os parâmetros \(\alpha\) e \(\lambda\) como \(2\). Inicialmente tomaremos \(n\) também igual a \(2\). Por fim, setamos o grid de valores que vamos trabalhar. Vale ressaltar que nossa distribuição tem suporte nos reais positivos.
## Setando valores
n <<- 2
alpha <<- 2
lambda <<- 2
s <- seq(1e-2, 10, by = 1e-2)Geramos, de forma pseudo-aleatório, \(10\) mil valores que seguem a distribuição WE para cada uma das n variáveis aleatórias e então somamos esses valores criando uma soma estocástica.
## Aproximação Empírica
# Geração de valores aleatórios pela soma de duas exponenciais
rwexp <- function(n,shape,scale)
{
V <- rexp(n, scale);
W <- rexp(n, (shape + 1)*scale);
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 de Váriáveis Aleatórias Contínuas
sn <- apply(x,1,sum)Finalmente, 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")
}
# Aproximação para a Densidade
par(mfrow = c(1,2), mar = c(3.2, 3.2, 1.6, 0.4), cex = 1.4)
plot(density(sn), type = 'l', lwd = 2, col = 1, lty = 1, xaxt = "n", yaxt = "n",
xlab = "", ylab = "", main = "", xlim=c(0,5), ylim=c(0,.65))
lines(s, densitysaddlepoint(s), lwd = 2, col = 1, lty = 2)
lines(s, dnorm(s,mean(sn),sd(sn)), lwd = 2, col = 1, lty = 3)
axis(2, seq(0, 0.65, length.out = 5), auxfunc(seq(0, 0.65, length.out = 5), 2))
axis(1, seq(0.01, 5, length.out = 5), auxfunc(seq(0.01, 5, length.out = 5), 2))
mtext("s", side = 1, line = 2.0, cex = 1.5)
mtext("Funções Densidade de Probabilidade", side = 2, line =2.0, cex = 1.5)
mtext("Aproximação Ponto de Sela (n=2)", side = 3, line =.4, cex = 1.5)
legend("topright", lty = 1:3, bty = "n", pch = c(20, 20, 20),
legend = c("Empírica", "Ponto de Sela","Normal"), cex = 1, inset = 0.02, lwd = 2)
# Aproximação para a Média
par(mar = c(3.2, 3.2, 1.6, 0.4), cex = 1.4)
plot(density(apply(x,1,mean)), type = 'l', lwd = 2, col = 1, lty = 1, xaxt = "n",
yaxt = "n", xlab = "", ylab = "", main = "", xlim=c(0,2.5), ylim=c(0,1.3))
lines(s/n, meansaddlepoint(s), lwd = 2, col = 1, lty = 2)
lines(s/n,dnorm(s/n,mean(apply(x,1,mean)),sd(apply(x,1,mean))),lwd=2,col=1,lty=3)
axis(2, seq(0, 1.3, length.out = 5), auxfunc(seq(0, 1.3, length.out = 5), 2))
axis(1, seq(0.01,2.5,length.out=5), auxfunc(seq(0.01,2.5,length.out=5),2))
mtext("s", side = 1, line = 2.0, cex = 1.5)
mtext("Funções Densidade da Média Amostral", side = 2, line =2.0, cex = 1.5)
mtext("Aproximação Ponto de Sela (n=2)", side = 3, line =.4, cex = 1.5)
legend("topright", lty = 1:3, bty = "n", pch = c(20, 20, 20),
legend = c("Empírica","Ponto de Sela","Normal"),cex=1,inset=0.02,lwd=2)
Figura 3: Aproximações Ponto de Sela Densidade (n = 2).
Analisando o painel esquerdo da Figura 3, é possível observar que a expansão ponto de sela para a densidade da soma estocástica tem um bom ajuste mesmo para um valor tão pequeno de \(n\). Observe que o fato da distribuição ser assimétrica e ter suporte positivo dificulta um ajuste adequado se utilizamos a distribuição normal. No painel direito também temos um bom ajuste da expansão ponto de sela para a densidade da média amostral.
# Aproximação para a Acumulada
par(mfrow = c(1,2), mar = c(3.2, 3.2, 1.6, 0.4), cex = 1.4)
plot(ecdf(sn), lwd = 2, col = 1, lty = 1, xaxt = "n", yaxt = "n",
xlab = "", ylab = "", main = "", xlim=c(0,5))
lines(s, danielssaddlepoint(s), lwd = 2, col = 1, lty = 2)
lines(s, LRsaddlepoint(s), lwd = 2, col = 1, lty = 3)
lines(s, pnorm(s,mean(sn),sd(sn)), lwd = 2, col = 1, lty = 4)
axis(2, seq(0, 1, length.out = 5), auxfunc(seq(0, 1, length.out = 5), 2))
axis(1, seq(0.01, 5, length.out = 5), auxfunc(seq(0.01, 5, length.out = 5), 2))
mtext("s", side = 1, line = 2.0, cex = 1.5)
mtext("Funções de Distribuição Acumulada", side = 2, line =2.0, cex = 1.5)
mtext("Aproximação Ponto de Sela (n=2)", side = 3, line =.4, cex = 1.5)
legend("bottomright", lty = 1:4, bty = "n", pch = c(20, 20, 20, 20),
legend = c("Empírica", "Daniels (1987)", "Lugannani e Rice (1980)", "Normal"),
cex = 1, inset = 0.02, lwd = 2)
# Comparação das Acumulada
par(mar = c(3.2, 3.2, 1.6, 0.4), cex = 1.4)
plot(s,danielssaddlepoint(s),type='l',lwd=2,col=1,lty=1,xaxt="n",yaxt="n",
xlab = "", ylab = "", main = "", xlim=c(0,5), ylim=c(0,1))
lines(s, LRsaddlepoint(s), lwd = 3, col = 1, lty = 3)
axis(2, seq(0, 1, length.out = 5), auxfunc(seq(0, 1, length.out = 5), 2))
axis(1, seq(0.01, 5, length.out = 5), auxfunc(seq(0.01, 5, length.out = 5), 2))
mtext("s", side = 1, line = 2.0, cex = 1.5)
mtext("Funções de Distribuição Acumulada", side = 2, line =2.0, cex = 1.5)
mtext("Comparação (n=2)", side = 3, line =.4, cex = 1.5)
legend("bottomright", lty = 2:3, bty = "n", pch = c(20, 20),
legend = c("Daniels (1987)", "Lugannani e Rice (1980)"),
cex = 1, inset = 0.02, lwd = 2)
Figura 4: Aproximações Ponto de Sela Distribuição Acumulada (n = 2).
No painel esquerdo da Figura 4, é possível verificar que tanto a distribuição acumulada de Daniels (1987) como a de Lugannani e Rice (1980) apresentam um ótimo ajuste, principalmente se levarmos em consideração o baixo valor de \(n\). No painel direito plotamos as duas acumuladas de Daniels (1987) e de Lugannani e Rice (1980) a fim de observar suas variabilidades. Fica claro que ambas as curvas estão sobrepostas, ficando indistinguíveis graficamente.
As Figuras 5 a 10 são construídas usando o código mostrado anteriormente, alterando apenas o valor de \(n\) para \(5\) e \(10\).
Figura 5: Aproximações Ponto de Sela para a Densidade de Probabilidade (n = 5).
Figura 6: Aproximações Ponto de Sela para a Densidade da Média Amostral (n = 5).
Figura 7: Aproximações Ponto de Sela para a Distribuição Acumulada (n = 5).
As Figuras 5, 6 e 7 mostram que as expansões ponto de sela estão descrevendo com excelente precisão a curva empírica. Em ambas, a distribuição normal começa a se aproximar da distribuição empírica.
Figura 8: Aproximações Ponto de Sela para a Densidade de Probabilidade (n = 10).
Figura 9: Aproximações Ponto de Sela para a Densidade da Média Amostral (n = 10).
Figura 10: Aproximações Ponto de Sela para a Distribuição Acumulada (n = 10).
Por fim, as Figuras 8, 9 e 10 mostram as curvas empíricas e as expansões ponto de sela praticamente iguais. Ao passo que o valor de \(n\) aumenta a distribuição normal se aproxima da distribuição empírica.
5 Conclusão
Este trabalho apresentou a distribuição exponencial ponderada criada por Gupta e Kundu (2009). Foi mostrado o comportamento da função densidade de probabilidade assim como o da função de risco. É notória a flexibilidade que o novo parâmetro de forma trouxe à distribuição base. Algumas propriedades foram apresentadas, entre elas, a função geradora de momentos, a função geratriz de momentos bem como seus quatro primeiros momentos que são respectivamente, média, variância, assimetria e curtose. Adiante foram apresentadas as expansões ponto de sela para a densidade da soma estocástica de variáveis aleatórias, para a densidade da média amostral, para a distribuição acumulada criada por Daniels (1987) e, por fim, para a distribuição acumulada criada por Lugannani e Rice (1980). Uma comparação foi feita a fim de avaliarmos quão boa é a qualidade das expansões ponto de sela se comparadas com a distribuição normal. Para \(n=2\) tivemos um bom ajuste proveniente das diversas expansões ponto de sela. Para \(n=5\) esse ajuste se tornou mais refinado. Já para \(n=10\) as expansões ponto de sela e a distribuição empírica estavam sobrepostas. Também foi possível observar que de acordo com que o valor de \(n\) aumentava, como esperado, a distribuição normal se aproximava cada vez mais da distribuição empírica. Porém, em todos os casos analisados é notória a superioridade das expansões ponto de sela para pequenos valores de \(n\). Vale ressaltar que tanto a distribuição de Daniels (1987) como a de Lugannani e Rice (1980), além de terem uma boa performance, seus valores são tão próximos a ponto de suas curvas se tornarem indistinguíveis. Dito isso, vale ressaltar que a distribuição criada por Lugannani e Rice (1980) é muito mais simples, logo fica como recomendação sua utilização em trabalhos futuros.
6 Referências
Azzalini, A. (1985). A Class of Distributions Which Includes the Normal Ones. Scandinavian Journal of Statistics, 12(2), 171-178.
Cordeiro, G. M. (1999). Introdução à teoria assintótica. Livro texto do 22º Colóquio Brasileiro de Matemática. IMPA.
Daniels, H. E. (1954). Saddlepoint approximations in statistics. Ann. Math.Statist. 25, 631–650.
Daniels, H. E. (1987). Tail probability approximations. Int. Stat. Rev., 55, 37-48.
Gupta, R. D. and Kundu, D. (2009). A new class of weighted exponential distributions. Statistics, 43:6, 621-634.
Lugannani, R. and Rice, S. (1980). Saddle point approximation for the distribution of the sum of independent random variables. Adv. Appl. Prob. 12, 475–490.
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.