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

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?

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$

$$\DeclareMathOperator{\E}{\mathbb{E}} M(t) = \E e^{t X}$$ $K(t)=\log M(t)$$\log$ $$e^{K(t)} = \int_{-\infty}^\infty e^{t x} f(x) \; dx =\int_{-\infty}^\infty \exp(tx+\log f(x) ) \; dx \\ = \int_{-\infty}^\infty \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_t$h(t,x) = -tx - \log f(x)$$h(t,x)$$x$$t$ $$h(t,x)=h(t,x_0) + h'(t,x_0)(x-x_0) +\frac12 h''(t,x_0) (x-x_0)^2 +\dotsm$$ $'$$x$ $$h'(t,x)=-t-\frac{\partial}{\partial x}\log f(x) \\ h''(t,x)= -\frac{\partial^2}{\partial x^2} \log f(x) > 0$$ $x_t$seja a solução para $h'(t,x_t)=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 $\dotsm$ , obtém $$e^{K(t)} \approx \int_{-\infty}^\infty \exp(-h(t,x_t)-\frac12 h''(t,x_t) (x-x_t)^2 ) \; dx \\ = e^{-h(t,x_t)} \int_{-\infty}^\infty e^{-\frac12 h''(t,x_t) (x-x_t)^2} \; dx$$ que é uma integral gaussiana, fornecendo $$e^{K(t)} \approx e^{-h(t,x_t)} \sqrt{\frac{2\pi}{h''(t,x_t)}}.$$ Isso fornece (uma primeira versão) da aproximação do ponto de sela como $$f(x_t) \approx \sqrt{\frac{h''(t,x_t)}{2\pi}} \exp(K(t) -t x_t) \\ \tag{*} \label{*}$$ 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 de$h'(t,x_t)=0$

$$t = -\frac{\partial}{\partial x_t} \log f(x_t).$$ $x_t$ $$\frac{\partial t}{\partial x_t} = -\frac{\partial^2}{\partial x_t^2} \log f(x_t) > 0$$ $t$$x_t$$x_t$$\frac{\partial}{\partial x_t} \log f(x_t)$ log f ( x t ) = K ( t ) - t x t - 1$\eqref{*}$ $$\log f(x_t) = K(t) -t x_t -\frac12 \log \frac{2\pi}{-\frac{\partial^2}{\partial x_t^2} \log f(x_t)}. \tag{**} \label{**}$$ xtxt∂logf(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. $x_t$$x_t$ $$\frac{\partial \log f(x_t)}{\partial x_t} \approx (K'(t)-x_t) \frac{\partial t}{\partial x_t} - t$$ $$0 \approx t + \frac{\partial \log f(x_t)}{\partial x_t} = (K'(t)-x_t) \frac{\partial t}{\partial x_t}$$ $t$$x_t$ $$K'(t) - x_t=0, \\ \tag{§} \label{§}$$

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 )$\eqref{*}$

$$h''(t,x_t) = -\frac{\partial^2 \log f(x_t)}{\partial x_t^2} \\ = -\frac{\partial}{\partial x_t} (\frac{\partial log f(x_t)}{\partial x_t} ) \\ = -\frac{\partial}{\partial x_t}(-t)= (\frac{\partial x_t}{\partial t})^{-1}$$ $K'(t)=x_t$ $$\frac{\partial x_t}{\partial t} = K''(t).$$ $$h''(t,x_t) = \frac1{K''(t)}$$ $f(x)$ $$f(x_t) \approx e^{K(t)- t x_t} \sqrt{\frac1{2\pi 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 .$x_t$$x_t$$t$

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 $n$$X_1, X_2, \dotsc, X_n$$n K(t)$

$$f(\bar{x}_t) = e^{nK(t) - n t \bar{x}_t} \sqrt{\frac{n}{2\pi 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)=\frac1{\sqrt{2\pi}} e^{-\frac12 x^2}$$ $M(t)=\exp(\frac12 t^2)$ $$K(t)=\frac12 t^2 \\ K'(t)=t \\ K''(t)=1$$ $t=x_t$ $$f(x_t) \approx e^{\frac12 t^2 -t x_t} \sqrt{\frac1{2\pi \cdot 1}} = \frac1{\sqrt{2\pi}} e^{-\frac12 x_t^2}$$

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): $X_1, X_2, \dotsc, X_n$$f$M ( t ) = 1

$$\hat{M}(t)= \frac1{n} \sum_{i=1}^n e^{t x_i}$$ K (t)=log M (t)log( H (T/n)n) K ˉ X (t)=nlog M (T/n)$\hat{K}(t) = \log \hat{M}(t)$$\log ( \hat{M}(t/n)^n )$ $$\hat{K}_{\bar{X}}(t) = n \log \hat{M}(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:

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.$\eqref{**}$

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!

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 ( λ ) = $x_1,\dots,x_n$$x\in R^d$ K '(λ)=y,yx1,...,x

$$\hat{K}(\lambda) = \frac{1}{n}\sum_{i=1}^{n}e^{\lambda^Tx_i},$$ $$\hat{K}'(\lambda) = y,$$ $y$$x_1,\dots,x_n$

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 .$\text{Gamma}(2, 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

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)

O ajuste é muito bom.

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 $\chi^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

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)