Como funciona a aproximação do ponto de sela?

38

Como funciona a aproximação do ponto de sela? Para que tipo de problema é bom?
(Sinta-se à vontade para usar um exemplo ou exemplos específicos a título ilustrativo)

Existem desvantagens, dificuldades, coisas a serem observadas ou armadilhas para os incautos?

Glen_b -Reinstate Monica
fonte

Respostas:

49

A aproximação do ponto de sela a uma função de densidade de probabilidade (funciona da mesma forma para funções de massa, mas falarei apenas aqui em termos de densidades) é uma aproximação surpreendentemente bem funcional, que pode ser vista como um refinamento no teorema do limite central. Portanto, ele só funcionará em ambientes onde existe um teorema do limite central, mas precisa de suposições mais fortes.

Começamos com a suposição de que a função geradora de momento existe e é duas vezes diferenciável. Isso implica, em particular, que todos os momentos existem. Seja uma variável aleatória com função de geração de momento (mgf) e cgf (função de geração cumulativa) (onde indica o logaritmo natural). No desenvolvimento, irei acompanhar de perto Ronald W. Butler: "Aproximações do ponto de sela com aplicativos" (CUP). Vamos desenvolver a aproximação do ponto de sela usando a aproximação de Laplace a uma determinada integral. Escrever X

M(t)=EetX
K(t)=registroM(t)registro
eK(t)=-etxf(x)dx=-exp(tx+registrof(x))dx=-exp(-h(t,x))dx
onde . Agora vamos Taylor expandir em considerando como uma constante. Isso dá onde ' denota diferenciação em relação a x . Observe que h '(t, x) = - t- \ frac {\ parcial} {\ parcial x} \ log f (x) \\ h' '(t, x) = - \ frac {\ parcial ^ 2} {\ parcial x ^ 2} \ log f (x)> 0 (a última desigualdade por suposição, pois é necessária para a aproximação funcionar). Deixe x_th(t,x)=-tx-registrof(x)h(t,x)xt
h(t,x)=h(t,x0 0)+h(t,x0 0)(x-x0 0)+1 12h(t,x0 0)(x-x0 0)2+
x
h(t,x)=-t-xregistrof(x)h(t,x)=-2x2registrof(x)>0 0
xtseja a solução para h(t,xt)=0 0 . Assumiremos que isso fornece um mínimo para h(t,x) em função de x . Usando essa expansão na integral e esquecendo a parte , obtém
eK(t)-exp(-h(t,xt)-1 12h(t,xt)(x-xt)2)dx=e-h(t,xt)-e-1 12h(t,xt)(x-xt)2dx
que é uma integral gaussiana, fornecendo
eK(t)e-h(t,xt)2πh(t,xt).
Isso fornece (uma primeira versão) da aproximação do ponto de sela como
(*)f(xt)h(t,xt)2πexp(K(t)-txt)
Observe que a aproximação tem a forma de uma família exponencial.

Agora precisamos fazer algum trabalho para obter isso de uma forma mais útil.

A partir de , obtemos Diferenciar isso em relação a fornece (pelas nossas suposições), de modo que a relação entre e é monótona, assim está bem definido. Precisamos de uma aproximação para . Para esse fim, resolvemos deh(t,xt)=0 0

t=-xtregistrof(xt).
xt
txt=-2xt2registrof(xt)>0 0
txtxtxtregistrof(xt) log f ( x t ) = K ( t ) - t x t - 1(*)
(**)registrof(xt)=K(t)-txt-1 12registro2π-2xt2registrof(xt).
xtxtlogf(xt) Assumindo que o último termo acima dependa apenas de , sua derivada em relação a é aproximadamente zero (voltaremos a comentar sobre isso), obtemos Até esta aproximação, temos então modo que e devam ser relacionados através da equação que é chamada de equação do ponto de sela. xtxt
registrof(xt)xt(K(t)-xt)txt-t
0 0t+registrof(xt)xt=(K(t)-xt)txt
txt
(§)K(t)-xt=0 0,

O que sentimos falta na determinação de é e que podemos encontrar por diferenciação implícita da equação do ponto de sela : O resultado é que (até nossa aproximação) Juntando tudo, temos a aproximação final do ponto de sela da densidade como H " ( t , x t ) = - 2 log de f ( x t )(*)

h(t,xt)=-2registrof(xt)xt2=-xt(euogf(xt)xt)=-xt(-t)=(xtt)-1 1
K(t)=xt
xtt=K(t).
h(t,xt)=1K(t)
f(x)
f(xt)eK(t)txt12πK(t).
xtxtt Agora, para usar isso praticamente, para aproximar a densidade em um ponto específico , resolvemos a equação do ponto de sela para que encontre .xtxtt

A aproximação saddlepoint é muitas vezes referido como uma aproximação para a densidade da média com base em observações iid . A função de geração cumulativa da média é simplesmente , portanto a aproximação do ponto de sela para a média se torna nX1,X2,,XnnK(t)

f(x¯t)=enK(t)ntx¯tn2πK(t)

Vejamos um primeiro exemplo. O que obtemos se tentarmos aproximar a densidade normal padrão O mgf é então para que a equação do ponto de sela seja e a aproximação do ponto de sela dê então, neste caso, a aproximação é exata.

f(x)=1 12πe-1 12x2
M(t)=exp(1 12t2)
K(t)=1 12t2K(t)=tK(t)=1 1
t=xt
f(xt)e1 12t2-txt1 12π1 1=1 12πe-1 12xt2

Vejamos um aplicativo muito diferente: Bootstrap no domínio de transformação, podemos executar bootstrap analiticamente usando a aproximação do ponto de sela à distribuição de bootstrap da média!

Suponha que temos iid distribuídos a partir de alguma densidade (no exemplo simulado, usaremos uma distribuição exponencial unitária). A partir da amostra, calculamos a função geradora de momento empírico e depois o cgf empírico . Precisamos do mgf empírico para a média que é e do cgf empírico para a média que usamos para construir uma aproximação do ponto de sela. No seguinte código R (versão R 3.2.3): X1 1,X2,...,XnfM ( t ) = 1

M^(t)=1 1nEu=1 1netxEu
K (t)=log M (t)log( H (T/n)n) K ˉ X (t)=nlog M (T/n)K^(t)=registroM^(t)registro(M^(t/n)n)
K^X¯(t)=nregistroM^(t/n)

set.seed(1234)
x  <-  rexp(10)

require(Deriv)   ### From CRAN
drule[["sexpmean"]]   <-  alist(t=sexpmean1(t))  # adding diff rules to 
                                                 # Deriv
drule[["sexpmean1"]]  <-  alist(t=sexpmean2(t))

###

make_ecgf_mean  <-   function(x)   {
    n  <-  length(x)
    sexpmean  <-  function(t) mean(exp(t*x))
    sexpmean1 <-  function(t) mean(x*exp(t*x))
    sexpmean2 <-  function(t) mean(x*x*exp(t*x))
    emgf  <-  function(t) sexpmean(t)
    ecgf  <-   function(t)  n * log( emgf(t/n) )
    ecgf1 <-   Deriv(ecgf)
    ecgf2 <-   Deriv(ecgf1)
    return( list(ecgf=Vectorize(ecgf),
                 ecgf1=Vectorize(ecgf1),
                 ecgf2 =Vectorize(ecgf2) )    )
}

### Now we need a function solving the saddlepoint equation and constructing
### the approximation:
###

make_spa <-  function(cumgenfun_list) {
    K  <- cumgenfun_list[[1]]
    K1 <- cumgenfun_list[[2]]
    K2 <- cumgenfun_list[[3]]
    # local function for solving the speq:
    solve_speq  <-  function(x) {
          # Returns saddle point!
          uniroot(function(s) K1(s)-x,lower=-100,
                  upper = 100, 
                  extendInt = "yes")$root
}
    # Function finding fhat for one specific x:
    fhat0  <- function(x) {
        # Solve saddlepoint equation:
        s  <-  solve_speq(x)
        # Calculating saddlepoint density value:
        (1/sqrt(2*pi*K2(s)))*exp(K(s)-s*x)
    }
    # Returning a vectorized version:
    return(Vectorize(fhat0))
} #end make_fhat

(Tentei escrever isso como um código geral que pode ser modificado facilmente para outros cgfs, mas o código ainda não é muito robusto ...)

Em seguida, usamos isso para uma amostra de dez observações independentes de uma distribuição exponencial unitária. Fazemos o bootstrap não paramétrico usual "à mão", plotamos o histograma de bootstrap resultante para a média e traçamos em excesso a aproximação do ponto de sela:

> ECGF  <- make_ecgf_mean(x)
> fhat  <-  make_spa(ECGF)
> fhat
function (x) 
{
    args <- lapply(as.list(match.call())[-1L], eval, parent.frame())
    names <- if (is.null(names(args))) 
        character(length(args))
    else names(args)
    dovec <- names %in% vectorize.args
    do.call("mapply", c(FUN = FUN, args[dovec], MoreArgs = list(args[!dovec]), 
        SIMPLIFY = SIMPLIFY, USE.NAMES = USE.NAMES))
}
<environment: 0x4e5a598>
> boots  <-  replicate(10000, mean(sample(x, length(x), replace=TRUE)), simplify=TRUE)
> boots  <-  replicate(10000, mean(sample(x, length(x), replace=TRUE)), simplify=TRUE)
> hist(boots, prob=TRUE)
> plot(fhat, from=0.001, to=2, col="red", add=TRUE)

Fornecendo a plotagem resultante:

aproximação do ponto de sela da distribuição de auto-inicialização

A aproximação parece ser bastante boa!

Poderíamos obter uma aproximação ainda melhor integrando a aproximação do ponto de sela e o redimensionamento:

> integrate(fhat, lower=0.1, upper=2)
1.026476 with absolute error < 9.7e-07

Agora, a função de distribuição cumulativa baseada nessa aproximação pode ser encontrada pela integração numérica, mas também é possível fazer uma aproximação direta do ponto de sela para isso. Mas isso é para outro post, isso é longo o suficiente.

Finalmente, alguns comentários foram deixados de fora do desenvolvimento acima. Em , fizemos uma aproximação essencialmente ignorando o terceiro termo. Por que podemos fazer isso? Uma observação é que, para a função de densidade normal, o termo deixado de fora não contribui em nada, de modo que a aproximação é exata. Portanto, como a aproximação do ponto de sela é um refinamento do teorema do limite central, estamos um pouco próximos do normal, portanto isso deve funcionar bem. Pode-se também olhar para exemplos específicos. Observando a aproximação do ponto de sela à distribuição de Poisson, observando o terceiro termo deixado de lado, neste caso, torna-se uma função trigamma, que de fato é bastante plana quando o argumento não é próximo de zero.(**)

Finalmente, por que o nome? O nome vem de uma derivação alternativa, usando técnicas de análise complexa. Mais tarde, podemos analisar isso, mas em outro post!

kjetil b halvorsen
fonte
4
O que você tem até agora é ótimo. O desenvolvimento lá é muito claro.
Glen_b -Replica Monica
11
kjetil Tentei consertar quatro pequenos erros de digitação 1. " No desenvolvimento, seguirei " 2. " necessário para a aproximação funcionar " 3. " O que sentimos falta agora " 4. " diferenciação implícita do sadlepoint ", mas ao fazê-lo parece que eu quebrei uma de suas equações - não faço idéia de como, pois não mudei nada além desses itens de texto (como você pode ver no histórico de edições). Eu o reverteria, mas como não consigo explicar como a correção desses erros causou um problema, não quero causar mais problemas. Me desculpe.
11
É possível que haja um bug mathJax ou um código no código de edição que leva a esse problema.
Glen_b -Reinstala Monica
11
@Christoph Hanck: Para obter uma aproximação em algum especifixo , você resolve a equação do ponto de sela para encontrar . (§) txt(§)t
Kjetil b halvorsen
2
Talvez valha a pena ressaltar que, quando o cgf empírico é usado, a aproximação resultante do ponto de sela é indefinida fora do casco convexo dos dados. Veja Feuerverger (1989) "Sobre a aproximação empírica do ponto de sela". Este deve ser o caso também no exemplo de inicialização acima.
Matteo Fasiolo
15

Aqui, eu expanda a resposta de kjetil e concentro-me nas situações em que a Função Geradora de Cumulantes (CGF) é desconhecida, mas pode ser estimada a partir dos dados , em que . O estimador CGF mais simples é provavelmente o de Davison e Hinkley (1988) que é o usado no exemplo de inicialização do kjetil. Esse estimador tem a desvantagem de que a equação resultante do ponto de sela só pode ser resolvida se , o ponto no qual queremos avaliar a densidade do ponto de sela, estiver dentro do casco convexo de . x R d K ( λ ) = 1x1 1,...,xnxRd K '(λ)=y,yx1,...,xn

K^(λ)=1 1nEu=1 1neλTxEu,
K^(λ)=y,
yx1 1,...,xn

Wong (1992) e Fasiolo et al. (2016) abordaram esse problema propondo dois estimadores de CGF alternativos, projetados de forma que a equação do ponto de sela possa ser resolvida para qualquer . A solução de Fasiolo et al. (2016), chamada de ESA de aproximação empírica empírica estendida, é implementada no pacote esaddle R e aqui vou dar alguns exemplos.y

Como um exemplo simples e univariado, considere usar o ESA para aproximar uma densidade .Gama(2,1 1)

library("devtools")
install_github("mfasiolo/esaddle")
library("esaddle")

########## Simulating data
x <- rgamma(1000, 2, 1)

# Fixing tuning parameter of ESA
decay <-  0.05

# Evaluating ESA at several point
xSeq <- seq(-2, 8, length.out = 200)
tmp <- dsaddle(y = xSeq, X = x, decay = decay, log = TRUE)

# Plotting true density, ESA and normal approximation
plot(xSeq, exp(tmp$llk), type = 'l', ylab = "Density", xlab = "x")
lines(xSeq, dgamma(xSeq, 2, 1), col = 3)
lines(xSeq, dnorm(xSeq, mean(x), sd(x)), col = 2)
suppressWarnings( rug(x) )
legend("topright", c("ESA", "Truth", "Gaussian"), col = c(1, 3, 2), lty = 1)

Este é o ajuste

insira a descrição da imagem aqui

Olhando para o tapete, fica claro que avaliamos a densidade da ESA fora do intervalo dos dados. Um exemplo mais desafiador é o seguinte gaussiano bivariado e distorcido.

# Function that evaluates the true density
dwarp <- function(x, alpha) {
  d <- length(alpha) + 1
  lik <- dnorm(x[ , 1], log = TRUE)
  tmp <- x[ , 1]^2
  for(ii in 2:d)
    lik <- lik + dnorm(x[ , ii] - alpha[ii-1]*tmp, log = TRUE)
  lik
}

# Function that simulates from true distribution
rwarp <- function(n = 1, alpha) {
  d <- length(alpha) + 1
  z <- matrix(rnorm(n*d), n, d)
  tmp <- z[ , 1]^2
  for(ii in 2:d) z[ , ii] <- z[ , ii] + alpha[ii-1]*tmp
  z
}

set.seed(64141)
# Creating 2d grid
m <- 50
expansion <- 1
x1 <- seq(-2, 3, length=m)* expansion; 
x2 <- seq(-3, 3, length=m) * expansion
x <- expand.grid(x1, x2) 

# Evaluating true density on grid
alpha <- 1
dw <- dwarp(x, alpha = alpha)

# Simulate random variables
X <- rwarp(1000, alpha = alpha)

# Evaluating ESA density
dwa <- dsaddle(as.matrix(x), X, decay = 0.1, log = FALSE)$llk

# Plotting true density
par(mfrow = c(1, 2))
plot(X, pch=".", col=1, ylim = c(min(x2), max(x2)), xlim = c(min(x1), max(x1)),
     main = "True density", xlab = expression(X[1]), ylab = expression(X[2]))
contour(x1, x2, matrix(dw, m, m), levels = quantile(as.vector(dw), seq(0.8, 0.995, length.out = 10)), col=2, add=T)

# Plotting ESA density
plot(X, pch=".",col=2, ylim = c(min(x2), max(x2)), xlim = c(min(x1), max(x1)),
     main = "ESA density", xlab = expression(X[1]), ylab = expression(X[2]))
contour(x1, x2, matrix(dwa, m, m), levels = quantile(as.vector(dwa), seq(0.8, 0.995, length.out = 10)), col=2, add=T)

insira a descrição da imagem aqui

O ajuste é muito bom.

Matteo Fasiolo
fonte
9

Graças à grande resposta de Kjetil, estou tentando criar um pequeno exemplo, que gostaria de discutir porque parece suscitar um ponto relevante:

Considere a . e seus derivados podem ser encontrados aqui e são reproduzidos nas funções no código abaixo. K ( t )χ2(m)K(t)

x <- seq(0.01,20,by=.1)
m <- 5

K  <- function(t,m) -1/2*m*log(1-2*t)
K1 <- function(t,m) m/(1-2*t)
K2 <- function(t,m) 2*m/(1-2*t)^2

saddlepointapproximation <- function(x) {
  t <- .5-m/(2*x)
  exp( K(t,m)-t*x )*sqrt( 1/(2*pi*K2(t,m)) )
}
plot( x, saddlepointapproximation(x), type="l", col="salmon", lwd=2)
lines(x, dchisq(x,df=m), col="lightgreen", lwd=2)

Isso produz

insira a descrição da imagem aqui

Obviamente, isso produz uma aproximação que acerta as características qualitativas da densidade, mas, como confirmado no comentário de Kjetil, não é uma densidade adequada, pois está acima da densidade exata em todos os lugares. O redimensionamento da aproximação da seguinte forma fornece o erro de aproximação quase insignificante plotado abaixo.

scalingconstant <- integrate(saddlepointapproximation, x[1], x[length(x)])$value

approximationerror_unscaled <- dchisq(x,df=m) - saddlepointapproximation(x)
approximationerror_scaled   <- dchisq(x,df=m) - saddlepointapproximation(x) /
                                                    scalingconstant

plot( x, approximationerror_unscaled, type="l", col="salmon", lwd=2)
lines(x, approximationerror_scaled,             col="blue",   lwd=2)

insira a descrição da imagem aqui

Christoph Hanck
fonte
11
É um recurso, a aproximação do ponto de sela não precisa ser integrada a um, mas costuma estar próxima. Pode ser redimensionado por integração numérica.
precisa saber é o seguinte
Poderia ser mais revelador traçar o erro relativo!
Kjetil b halvorsen
approximationerror_unscaled/approximationerror_scaledAcontece que a pairar em torno 25,90798
Christoph Hanck