Jak działa przybliżanie saddlepoint?

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)

M (t) = E e^{t X}

$\DeclareMathOperator{\E}{\mathbb{E}} M(t) = \E e^{t X}$ i cgf (funkcja generowania skumulowanego)

K (t) = \log M (t)

$K(t)=\log M(t)$ (gdzie

\log

$\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) d x = \int_{- \infty}^{\infty} \exp (t x + \log f (x)) d x = \int_{- \infty}^{\infty} \exp (- h (t, x)) d x

$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) = - t x - \log f (x)

$h(t,x) = -tx - \log f(x)$ . Teraz Taylor rozszerzy

h (t, x)

$h(t,x)$ w

x

$x$ traktując

t

$t$ jako stałą. Daje to

h (t, x) = h (t, x_{0}) + h^{'} (t, x_{0}) (x - x_{0}) + \frac{1}{2} h^{″} (t, x_{0}) (x - x_{0})^{2} + \dots

$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

$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

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

$x_t$ będzie rozwiązaniem dla

h^{'} (t, x_{t}) = 0

$h'(t,x_t)=0$ . Zakładamy, że daje to minimum dla

h (t, x)

$h(t,x)$ w funkcji

x

$x$ . Używając tego rozszerzenia w całce i zapominając oczęści

\dots

$\dotsm$ , daje

e^{K (t)} \approx \int_{- \infty}^{\infty} \exp (- h (t, x_{t}) - \frac{1}{2} h^{″} (t, x_{t}) (x - x_{t})^{2}) d x = e^{- h (t, x_{t})} \int_{- \infty}^{\infty} e^{- \frac{1}{2} h^{″} (t, x_{t}) (x - x_{t})^{2}} d x

$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 π}{h^{″} (t, x_{t})}} .

$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

\begin{matrix} (*) & f (x_{t}) \approx \sqrt{\frac{h^{″} (t, x_{t})}{2 π}} \exp (K (t) - t x_{t}) \end{matrix}

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

$t = -\frac{\partial}{\partial x_t} \log f(x_t).$ Zróżnicowanie tego względem

x_{t}

$x_t$ daje

\frac{\partial t}{\partial x_{t}} = - \frac{\partial^{2}}{\partial x_{t}^{2}} \log f (x_{t}) > 0

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

$t$ i

x_{t}

$x_t$ znaczy monotoniczne tak

x_{t}

$x_t$ jest dobrze zdefiniowana. Potrzebujemy przybliżenia do

\frac{\partial}{\partial x_{t}} \log f (x_{t})

$\frac{\partial}{\partial x_t} \log f(x_t)$ . W tym celu otrzymujemy rozwiązanie z

(*)

$\eqref{*}$

\begin{matrix} (**) & \log f (x_{t}) = K (t) - t x_{t} - \frac{1}{2} \log \frac{2 π}{- \frac{\partial^{2}}{\partial x_{t}^{2}} \log f (x_{t})} . \end{matrix}

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

$x_t$ , więc jego pochodna względem

x_{t}

$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

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

$0 \approx t + \frac{\partial \log f(x_t)}{\partial x_t} = (K'(t)-x_t) \frac{\partial t}{\partial x_t}$ , aby

t

$t$ i

x_{t}

$x_t$ były powiązane przez równanie

\begin{matrix} (§) & K^{'} (t) - x_{t} = 0, \end{matrix}

$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 l o g f (x_{t})}{\partial x_{t}}) = - \frac{\partial}{\partial x_{t}} (- t) = (\frac{\partial x_{t}}{\partial t})^{- 1}

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

$K'(t)=x_t$

\frac{\partial x_{t}}{\partial t} = K^{″} (t) .

$\frac{\partial x_t}{\partial t} = K''(t).$

h^{″} (t, x_{t}) = \frac{1}{K^{″} (t)}

$h''(t,x_t) = \frac1{K''(t)}$

f (x)

$f(x)$

f (x_{t}) \approx e^{K (t) - t x_{t}} \sqrt{\frac{1}{2 π K^{″} (t)}} .

$f(x_t) \approx e^{K(t)- t x_t} \sqrt{\frac1{2\pi K''(t)}}.$

x_{t}

$x_t$

x_{t}

$x_t$

t

$t$

$n$ $X_1, X_2, \dotsc, X_n$ $n K(t)$

f ({\bar{x}}_{t}) = e^{n K (t) - n t {\bar{x}}_{t}} \sqrt{\frac{n}{2 π K^{″} (t)}}

$f(\bar{x}_t) = e^{nK(t) - n t \bar{x}_t} \sqrt{\frac{n}{2\pi K''(t)}}$

f (x) = \frac{1}{\sqrt{2 π}} e^{- \frac{1}{2} x^{2}}

$f(x)=\frac1{\sqrt{2\pi}} e^{-\frac12 x^2}$

M (t) = \exp (\frac{1}{2} t^{2})

$M(t)=\exp(\frac12 t^2)$

K (t) = \frac{1}{2} t^{2} K^{'} (t) = t K^{″} (t) = 1

$K(t)=\frac12 t^2 \\ K'(t)=t \\ K''(t)=1$

t = x_{t}

$t=x_t$

f (x_{t}) \approx e^{\frac{1}{2} t^{2} - t x_{t}} \sqrt{\frac{1}{2 π \cdot 1}} = \frac{1}{\sqrt{2 π}} e^{- \frac{1}{2} x_{t}^{2}}

$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) = \frac{1}{n} \sum_{i = 1}^{n} e^{t x_{i}}

$\hat{M}(t)= \frac1{n} \sum_{i=1}^n e^{t x_i}$

\hat{K} (t) = \log \hat{M} (t)

$\hat{K}(t) = \log \hat{M}(t)$

\log (\hat{M} (t / n)^{n})

$\log ( \hat{M}(t/n)^n )$

{\hat{K}}_{\bar{X}} (t) = n \log \hat{M} (t / 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!

kjetil b halvorsen
źródło

To, co do tej pory masz, jest świetne. Rozwój tam jest bardzo wyraźny.

Glen_b

kjetil Próbowałem naprawić cztery małe literówki 1. „ W rozwoju będę śledzić ” 2. „ Potrzebne do przybliżenia do działania ” 3. „ To, czego teraz brakuje ” 4. „ Ukryte różnicowanie punktu siodłowego ”, ale robiąc to wygląda na to, że złamałem jedno z twoich równań - nie mam pojęcia, jak to zrobić, ponieważ nie zmieniłem nic oprócz tych elementów tekstowych (jak widać z historii edycji). Wycofam go, ale ponieważ nie potrafię wyjaśnić, w jaki sposób naprawienie tych błędów spowodowało problem, nie chcę powodować dalszych problemów. Przepraszam. (Wyglądało na to, że się zepsuł, gdy tylko otworzyłem sesję edycji)

Glen_b

Możliwe, że w kodzie edycji występuje błąd mathJax lub błąd, który prowadzi do tego problemu.

Glen_b

x_{t}

$x_t$

(§)

$\eqref{§}$

t

$t$

Być może warto zauważyć, że gdy stosuje się empiryczne cgf, wynikowe przybliżenie punktu saddlepoint jest niezdefiniowane poza wypukłym kadłubem danych. Zobacz Feuerverger (1989) „O empirycznym przybliżeniu punktu siodłowego”. Tak powinno być również w powyższym przykładzie bootstrap.

Matteo Fasiolo

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

\hat{K} (λ) = \frac{1}{n} \sum_{i = 1}^{n} e^{λ^{T} x_{i}},

$\hat{K}(\lambda) = \frac{1}{n}\sum_{i=1}^{n}e^{\lambda^Tx_i},$

{\hat{K}}^{'} (λ) = y,

$\hat{K}'(\lambda) = y,$

y

$y$

x_{1}, \dots, x_{n}

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

Matteo Fasiolo
źródło

Jak działa przybliżanie saddlepoint?

Odpowiedzi: