Expansões(Andre) Ponto Sela

Prévia do material em texto

Expansão Ponto de Sela
Vamos agora aproximar a soma de variáveis aleatórias com certa distribuição através da metodologia ponto
de sela. Ao contrário da expansão de Edgeworth, nesta é necessário que a distribuição possua função geratriz
de cumulantes (f.g.c.). Com isso, não podemos usar a expansão ponto de sela na distribuição Weibull. Vamos
usar agora a distribuição Gamma.
A expansão ponto de sela para aproximar a função de densidade de Sn é dada por:
fsn(s) =
exp(nφ(λ̂)− sλ̂)√
2nπφ(2)(λ̂)
[
1 +
(
3ρ4(λ̂)− 5ρ3(λ̂)2
24n
)
+On−2
]
, (5)
em que φ(λ) é a função geratriz de cumulantes; φ(r)(λ̂) = d(r)φ(λ)/dλr; ρr(λ) = φ(r)(λ)/φ(2)(λ)r/2; e λ̂ é a
estimativa de máxima-verossimilhança que satisfaz φ(1)(λ) = s/n.
Distribuição Gamma
A distribuição Gamma, assim, como a Weibull possuem dois parâmetros. Sua função de densidade (fdp) é:
f(x;α, β) = αΓ(β) (αx)
β−1e−αx, com α, β > 0. (6)
# A exponencial é tambem um caso especial da gamma
dgamma(7,shape = 1, rate = 3) == dexp(7, rate = 3)
## [1] TRUE
# plot gammapara alguns parametros
curve(dgamma(x, shape = 3, rate=7), xlim = c(0,7), lwd = 2,
main = "Distribuição Gamma para alguns parâmetros",
xlab = "x", ylab = "Densidade")
curve(dgamma(x, shape = 3, rate=5), xlim = c(0,7), lwd = 2, add = T)
curve(dgamma(x, shape = 4, rate=2), xlim = c(0,7), lwd = 2, add = T)
curve(dgamma(x, shape = 6, rate=2), xlim = c(0,7), lwd = 2, add = T)
text(1.15,1.6, expression(paste(alpha==7," ", "e", " ", beta==3)))
text(1.37,1, expression(paste(alpha==5," ", "e", " ", beta==3)))
text(2,.5, expression(paste(alpha==2," ", "e", " ", beta==4)))
text(4,.34, expression(paste(alpha==2," ", "e", " ", beta==6)))
8
0 1 2 3 4 5 6 7
0.
0
0.
5
1.
0
1.
5
Distribuição Gamma para alguns parâmetros
x
D
en
si
da
de
α = 7 e β = 3
α = 5 e β = 3
α = 2 e β = 4
α = 2 e β = 6
Fica visível no gráfico que α e β na equação (6) são, respectivamente, os parâmetros de escala e de forma.
A função geradora de momentos (fgm) da fdp (6) é dada por:
M(λ) =
(
1− λ
α
)−β
. (7)
Como a fgc é o logaritmo da fgm, temos:
φ(λ) = log
((
1− λ
α
)−β)
=− β log
(
1− λ
α
)
(8)
A r-ésima derivada de (8), φ(r)(λ), produz o r-ésimo cumulante. Assim,
φ(1)(λ) = β
α− λ
, (9)
φ(2)(λ) = β(α− λ)2 , (10)
φ(3)(λ) = 2β(α− λ)3 , (11)
φ(4)(λ) = 6β(α− λ)4 . (12)
Fazendo, φ(1)(λ) = s/n, na equação (9)
β
α− λ
= s
n
.
9
Tendo como solução para λ
λ̂ = α− βn
s
. (13)
Desta maneira, λ nas equações (8), (10), (11) e (12), é substituído por λ̂ da equação (13).
# Expansao ponto de sela
ponto_sela <- function(s, n, lambda, K_lambda, k2_lambda, rho3, rho4){
M.l <- (3 * rho4 - 5*rho3^2) / (24 * n )
qa <- exp(n*K_lambda - s*lambda)
qs <- sqrt(2*n*pi*k2_lambda)
f_sn <- (qa / qs) * (1 + M.l)
f_sn
}
lambda <- function(alpha, beta, n, s){
qa <- alpha - (beta*n)/s
qa
}
# Funcao fgc no ponto em que o primeiro cumlante eh igual a s/n
K_lambda <- function(alpha, beta, lambda){
qa <- -beta * log(1 - lambda/alpha)
}
# segundo cumulante no ponto lambda = s/n
K2_lambda <- function(alpha, beta, lambda){
qa <- beta/ (alpha - lambda)^2
qa
}
# terceiro cumulante no ponto lambda = s/n
K3_lambda <- function(alpha, beta, lambda){
qa <- 2*beta/ (alpha - lambda)^3
qa
}
# quarto cumulante no ponto lambda = s/n
K4_lambda <- function(alpha, beta, lambda){
qa <- 6*beta/ (alpha - lambda)^4
qa
}
s.gama <- seq(from = 0.2, by=.1, to=40)
beta.gama <- 4
alfa.gama <- 2
# obtendo lambda para n = 5
la5 <- lambda(alpha = alfa.gama, beta = beta.gama, n=5, s=s.gama)
# obtendo o valor da fgc no em lambda
10
k.l5 <- K_lambda(alpha = alfa.gama, beta = beta.gama, lambda = la5)
# obtendo o segundo cumulante
k2.l5 <- K2_lambda(alpha = alfa.gama, beta = beta.gama, lambda = la5)
# obtendo o terceiro cumulante
k3.l5 <- K3_lambda(alpha = alfa.gama, beta = beta.gama, lambda = la5)
# obtendo o quarto cumulante
k4.l5 <- K4_lambda(alpha = alfa.gama, beta = beta.gama, lambda = la5)
# padronizando os cumulantes
rho3.l5 <- k3.l5 / k2.l5^(3/2)
rho4.l5 <- k4.l5 / k2.l5^(4/2)
# obtendo lambda para n = 10
la10 <- lambda(alpha = alfa.gama, beta = beta.gama, n=10, s=s.gama)
# obtendo o valor da fgc no em lambda
k.l10 <- K_lambda(alpha = alfa.gama, beta = beta.gama, lambda = la10)
# obtendo o segundo cumulante
k2.l10 <- K2_lambda(alpha = alfa.gama, beta = beta.gama, lambda = la10)
# obtendo o terceiro cumulante
k3.l10 <- K3_lambda(alpha = alfa.gama, beta = beta.gama, lambda = la10)
# obtendo o quarto cumulante
k4.l10 <- K4_lambda(alpha = alfa.gama, beta = beta.gama, lambda = la10)
# padronizando os cumulantes
rho3.l10 <- k3.l10 / k2.l10^(3/2)
rho4.l10 <- k4.l10 / k2.l10^(4/2)
# ponto de sela
point.sela5 <- ponto_sela(s=s.gama,n=5,lambda = la5, K_lambda = k.l5,
k2_lambda = k2.l5, rho3 = rho3.l5, rho4 = rho4.l5 )
point.sela10 <- ponto_sela(s=s.gama,n=10,lambda = la10, K_lambda = k.l10,
k2_lambda = k2.l10, rho3 = rho3.l10, rho4 = rho4.l10 )
## Gerando numeros aleatorios gama
# n = 5
sngama5 <- matrix(rgamma(200000, shape = beta.gama, rate = alfa.gama),
byrow = T, ncol = 5 )
# n = 10
11
sngama10 <- matrix(rgamma(400000, shape = beta.gama, rate = alfa.gama),
byrow = T, ncol = 10 )
soma_gama5 <- rowSums(sngama5)
soma_gama10 <- rowSums(sngama10)
plot(density(soma_gama5), type = "l", lwd = 2, main = expression(
paste("Aproximação Ponto de Sela de"," ", S[n], " ", "(n=5)")) ,
ylab = "Densidade", xlab = expression(paste(S[n] == Y[1]+"..."+Y[n])))
lines(s.gama, point.sela5, col = 2, lwd = 2)
legend("topright", legend = c("Empírica", "Ponto de Sela"),
col = c(1,2), lwd = 2 , bty ="n")
5 10 15 20
0.
00
0.
05
0.
10
0.
15
Aproximação Ponto de Sela de Sn (n=5)
Sn = Y1 + ... + Yn
D
en
si
da
de
Empírica
Ponto de Sela
plot(density(soma_gama10), type = "l", lwd = 2, main = expression(
paste("Aproximação Ponto de Sela de"," ", S[n], " ", "(n=10)")) ,
ylab = "Densidade", xlab = expression(paste(S[n] == Y[1]+"..."+Y[n])), ylim=c(0,.13) )
lines(s.gama, point.sela10, col = 2, lwd = 2)
legend("topright", legend = c("Empírica", "Ponto de Sela"),
col = c(1,2), lwd = 2 , bty ="n")
12
10 15 20 25 30 35
0.
00
0.
04
0.
08
0.
12
Aproximação Ponto de Sela de Sn (n=10)
Sn = Y1 + ... + Yn
D
en
si
da
de
Empírica
Ponto de Sela
13
	Expansão Ponto de Sela
	Distribuição Gamma