Jak działa przybliżanie saddlepoint?

Jak działa przybliżanie saddlepoint? Dla jakiego rodzaju problemu jest to dobre?
(Możesz użyć konkretnego przykładu lub przykładów jako ilustracji)

Czy są jakieś wady, trudności, rzeczy, na które należy uważać, lub pułapki na nieostrożnych?

Przybliżenie punktu sadd do funkcji gęstości prawdopodobieństwa (działa podobnie w przypadku funkcji masy, ale będę mówić tutaj tylko w odniesieniu do gęstości) jest zaskakująco dobrze działającym przybliżeniem, które można postrzegać jako udoskonalenie centralnego twierdzenia o granicy. Będzie więc działał tylko w ustawieniach, w których istnieje centralne twierdzenie o limicie, ale wymaga silniejszych założeń.

Zaczynamy od założenia, że ​​funkcja generowania momentu istnieje i jest dwukrotnie rozróżnialna. Oznacza to w szczególności, że wszystkie chwile istnieją. Niech $X$ będzie zmienną losową z funkcją generowania momentu (mgf)

$$\DeclareMathOperator{\E}{\mathbb{E}} M(t) = \E e^{t X}$$ i cgf (funkcja generowania skumulowanego) $K(t)=\log M(t)$ (gdzie $\log$oznacza logarytm naturalny). W trakcie rozwoju będę uważnie śledzić Ronalda W Butlera: „Saddlepoint Approximations with Applications” (CUP). Rozwiniemy przybliżenie punktu siodłowego za pomocą przybliżenia Laplace'a do pewnej całki. Napisz $$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$$ gdzie $h(t,x) = -tx - \log f(x)$ . Teraz Taylor rozszerzy$h(t,x)$ w$x$ traktując$t$ jako stałą. Daje to $$h(t,x)=h(t,x_0) + h'(t,x_0)(x-x_0) +\frac12 h''(t,x_0) (x-x_0)^2 +\dotsm$$ gdzie$'$ Oznacza różnicowanie w odniesieniu do$x$ . Zauważ, że $$h'(t,x)=-t-\frac{\partial}{\partial x}\log f(x) \\ h''(t,x)= -\frac{\partial^2}{\partial x^2} \log f(x) > 0$$ (ostatnia nierówność z założenia, ponieważ jest to potrzebne do przybliżenia, aby zadziałało). Niech$x_t$będzie rozwiązaniem dla$h'(t,x_t)=0$. Zakładamy, że daje to minimum dla $h(t,x)$w funkcji$x$. Używając tego rozszerzenia w całce i zapominając oczęści$\dotsm$, daje $$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$$ która jest całką Gaussa, dającą $$e^{K(t)} \approx e^{-h(t,x_t)} \sqrt{\frac{2\pi}{h''(t,x_t)}}.$$ Daje to (pierwsza wersja) przybliżenie punktu siodłowego jako $$f(x_t) \approx \sqrt{\frac{h''(t,x_t)}{2\pi}} \exp(K(t) -t x_t) \\ \tag{*} \label{*}$$ Należy zauważyć, że przybliżenie ma postać wykładniczej rodziny.

Teraz musimy trochę popracować, aby uzyskać to w bardziej użytecznej formie.

Z $h'(t,x_t)=0$ otrzymujemy

$$t = -\frac{\partial}{\partial x_t} \log f(x_t).$$ Zróżnicowanie tego względem$x_t$daje $$\frac{\partial t}{\partial x_t} = -\frac{\partial^2}{\partial x_t^2} \log f(x_t) > 0$$ (o przyjętych założeń), a więc stosunek między$t$i$x_t$znaczy monotoniczne tak$x_t$jest dobrze zdefiniowana. Potrzebujemy przybliżenia do$\frac{\partial}{\partial x_t} \log f(x_t)$. W tym celu otrzymujemy rozwiązanie z 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{**}$$ Zakładając, że powyższy ostatni termin słabo zależy tylko od$x_t$, więc jego pochodna względem$x_t$wynosi w przybliżeniu zero (wrócimy do komentowania tego), otrzymujemy $$\frac{\partial \log f(x_t)}{\partial x_t} \approx (K'(t)-x_t) \frac{\partial t}{\partial x_t} - t$$ Do tego przybliżenia mamy wtedy $$0 \approx t + \frac{\partial \log f(x_t)}{\partial x_t} = (K'(t)-x_t) \frac{\partial t}{\partial x_t}$$ , aby$t$i$x_t$były powiązane przez równanie $$K'(t) - x_t=0, \\ \tag{§} \label{§}$$ który nazywa się równaniem saddlepoint.

$\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)}}.$$ $x_t$$x_t$$t$

$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)}}$$

$$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}$$

Spójrzmy na zupełnie inną aplikację: Bootstrap w domenie transformacji, możemy wykonać bootstrap analitycznie, używając przybliżenia saddlepoint do rozkładu bootstrap średniej!

$X_1, X_2, \dotsc, X_n$$f$

$$\hat{M}(t)= \frac1{n} \sum_{i=1}^n e^{t x_i}$$ $\hat{K}(t) = \log \hat{M}(t)$$\log ( \hat{M}(t/n)^n )$ $$\hat{K}_{\bar{X}}(t) = n \log \hat{M}(t/n)$$ których używamy do konstruowania przybliżenia punktu saddlepoint. Poniżej niektóre kod R (wersja R 3.2.3):

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

(Próbowałem napisać ten kod jako ogólny kod, który można łatwo modyfikować dla innych programów cgfs, ale kod nadal nie jest zbyt solidny ...)

Następnie używamy tego do próby dziesięciu niezależnych obserwacji z jednostkowego rozkładu wykładniczego. Wykonujemy zwykłe nieparametryczne ładowanie początkowe „ręcznie”, wykreślamy wynikowy histogram ładowania dla średniej i zastępujemy przybliżenie punktu Saddlepoint:

> 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)

Dając wynikowy wykres:

Przybliżenie wydaje się raczej dobre!

Możemy uzyskać jeszcze lepsze przybliżenie, integrując przybliżenie punktu saddpap i przeskalowanie:

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

Teraz funkcję rozkładu skumulowanego opartą na tym przybliżeniu można znaleźć za pomocą całkowania numerycznego, ale możliwe jest również bezpośrednie przybliżenie tego punktu. Ale to jest na inny post, to wystarczy.

$\eqref{**}$

Wreszcie, dlaczego nazwa? Nazwa pochodzi od alternatywnej pochodnej, wykorzystującej techniki analizy złożonej. Później możemy przyjrzeć się temu, ale w innym poście!

$x_1,\dots,x_n$$x\in R^d$

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

$y$

$\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)

To pasuje

Patrząc na dywan, widać, że oceniliśmy gęstość ESA poza zakresem danych. Bardziej wymagającym przykładem jest następujący wypaczony dwuwymiarowy gaussowski.

# 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)

Dopasowanie jest całkiem dobre.

Dzięki świetnej odpowiedzi Kjetil sam próbuję wymyślić mały przykład, który chciałbym omówić, ponieważ wydaje się, że porusza on istotną kwestię:

$\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)

To produkuje

To oczywiście daje przybliżenie, które poprawia cechy jakościowe gęstości, ale, jak potwierdzono w komentarzu Kjetila, nie jest właściwą gęstością, ponieważ jest wszędzie powyżej dokładnej gęstości. Ponowne skalowanie aproksymacji w następujący sposób daje niemal pomijalny błąd aproksymacji przedstawiony poniżej.

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)