Wie funktioniert die Sattelpunktnäherung?

38

Wie funktioniert die Sattelpunktnäherung? Für was für ein Problem ist es gut?
(Fühlen Sie sich frei, ein bestimmtes Beispiel oder Beispiele zur Veranschaulichung zu verwenden)

Gibt es Nachteile, Schwierigkeiten, Dinge, auf die man achten muss, oder Fallen für die Unachtsamen?

Glen_b
quelle

Antworten:

48

Die Sattelpunktnäherung an eine Wahrscheinlichkeitsdichtefunktion (sie funktioniert ebenfalls für Massenfunktionen, aber ich werde hier nur auf die Dichten eingehen) ist eine überraschend gut funktionierende Näherung, die als Verfeinerung des zentralen Grenzwertsatzes angesehen werden kann. Es funktioniert also nur in Umgebungen, in denen es einen zentralen Grenzwertsatz gibt, es sind jedoch strengere Annahmen erforderlich.

Wir gehen davon aus, dass die momenterzeugende Funktion existiert und doppelt differenzierbar ist. Dies impliziert insbesondere, dass alle Momente existieren. Sei eine Zufallsvariable mit Momenterzeugungsfunktion (mgf) und cgf (kumulative Erzeugungsfunktion) (wobei den natürlichen Logarithmus bezeichnet). In der Entwicklung werde ich Ronald W. Butler genau verfolgen: "Sattelpunktapproximationen mit Anwendungen" (CUP). Wir werden die Sattelpunktnäherung unter Verwendung der Laplace-Näherung auf ein bestimmtes Integral entwickeln. Schreiben X

M(t)=EetX
K(t)=logM(t)log
eK(t)=etxf(x)dx=exp(tx+logf(x))dx=exp(h(t,x))dx
wobei . Nun wollen wir Taylor in , wobei als Konstante betrachtet wird. Dies ergibt wobei ' bezeichnet die Differenzierung bezüglich x . Man beachte, dass h '(t, x) = - t - \ frac {\ partial} {\ partial x} \ log f (x) \\ h' '(t, x) = - \ frac {\ partial ^ 2} {\ partial x ^ 2} \ log f (x)> 0 (die letzte Ungleichung nach Annahme, wie sie für die Annäherung benötigt wird). Sei x_th(t,x)=txlogf(x)h(t,x)xt
h(t,x)=h(t,x0)+h(t,x0)(xx0)+12h(t,x0)(xx0)2+
x
h(t,x)=txlogf(x)h(t,x)=2x2logf(x)>0
xtsei die Lösung für h(t,xt)=0 . Wir nehmen an, dass dies ein Minimum für h(t,x) als Funktion von x . Wenn Sie diese Erweiterung im Integral verwenden und den Teil vergessen , erhalten Sie
eK(t)exp(h(t,xt)12h(t,xt)(xxt)2)dx=eh(t,xt)e12h(t,xt)(xxt)2dx
ist ein Gaußsches Integral und ergibt
eK(t)eh(t,xt)2πh(t,xt).
Dies ergibt (eine erste Version) der Sattelpunktnäherung als
(*)f(xt)h(t,xt)2πexp(K(t)txt)
Beachten Sie, dass die Approximation die Form einer Exponentialfamilie hat.

Jetzt müssen wir einige Arbeiten erledigen, um dies in eine nützlichere Form zu bringen.

Aus wir Wenn man dies in Bezug auf differenziert, erhält man (nach unseren Annahmen), Die Beziehung zwischen und ist also monoton, also ist gut definiert. Wir brauchen eine Annäherung an . Zu diesem Zweck lösen wir aush(t,xt)=0

t=xtlogf(xt).
xt
txt=2xt2logf(xt)>0
txtxtxtlogf(xt)(*)
(**)logf(xt)=K(t)txt12log2π2xt2logf(xt).
Unter der Annahme, dass der letzte Term oben nur schwach von abhängt , so dass seine Ableitung in Bezug auf ungefähr Null ist (wir werden darauf zurückkommen, um dies zu kommentieren), erhalten wir Bis zu dieser Annäherung haben wir dann diese so dass und durch die Gleichung in Beziehung gesetzt werden müssen , die Sattelpunktgleichung genannt wird. xtxt
logf(xt)xt(K(t)xt)txtt
0t+logf(xt)xt=(K(t)xt)txt
txt
(§)K(t)xt=0,

Was wir jetzt bei der Bestimmung von vermissen, ist und das können wir durch implizite Differenzierung der Sattelpunktgleichung finden : Das Ergebnis ist, dass (bis zu unserer Näherung) Wenn wir alles zusammenfassen, haben wir die endgültige Sattelpunktnäherung der Dichte als (*)

h(t,xt)=2logf(xt)xt2=xt(logf(xt)xt)=xt(t)=(xtt)1
K(t)=xt
xtt=K(t).
h(t,xt)=1K(t)
f(x)
f(xt)eK(t)txt12πK(t).
Um dies praktisch anzuwenden, um die Dichte an einem bestimmten Punkt anzunähern, lösen wir die Sattelpunktgleichung für diesen , um zu finden .xtxtt

Die Sattelpunktapproximation wird oft als eine Annäherung an die Dichte des Mittelwerts angegeben , basierend auf IId Beobachtungen . Die kumulative Erzeugungsfunktion des Mittelwerts ist einfach , so dass die Sattelpunktnäherung für den Mittelwert zu nX1,X2,,XnnK(t)

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

Schauen wir uns ein erstes Beispiel an. Was erhalten wir, wenn wir versuchen, die normale Standarddichte anzunähern. Die mgf ist also so dass die Sattelpunktgleichung und die Sattelpunktnäherung ergibt also in diesem Fall ist die Annäherung genau.

f(x)=12πe12x2
M(t)=exp(12t2)
K(t)=12t2K(t)=tK(t)=1
t=xt
f(xt)e12t2txt12π1=12πe12xt2

Betrachten wir eine ganz andere Anwendung: Bootstrap in der Transformationsdomäne, wir können Bootstrapping analytisch durchführen, indem wir die Sattelpunktannäherung an die Bootstrap-Verteilung des Mittelwerts verwenden!

Angenommen, wir haben iid von einer Dichte (im simulierten Beispiel verwenden wir eine Exponentialverteilung in Einheiten). Aus der Stichprobe berechnen wir die empirische Momenterzeugungsfunktion und dann die empirische cgf . Wir benötigen die empirische mgf für den Mittelwert, der und die empirische cgf für den Mittelwert dem wir eine Sattelpunktnäherung konstruieren. Im Folgenden einige R-Code (R-Version 3.2.3): X1,X2,,Xnf

M^(t)=1ni=1netxi
K^(t)=logM^(t)log(M^(t/n)n)
K^X¯(t)=nlogM^(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

(Ich habe versucht, dies als allgemeinen Code zu schreiben, der leicht für andere CGFS geändert werden kann, aber der Code ist immer noch nicht sehr robust ...)

Dann verwenden wir dies für eine Stichprobe von zehn unabhängigen Beobachtungen aus einer Exponentialverteilung. Wir machen das übliche nichtparametrische Bootstrapping "von Hand", zeichnen das resultierende Bootstrap-Histogramm für den Mittelwert und überzeichnen die Sattelpunkt-Approximation:

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

Geben Sie die resultierende Handlung:

Sattelpunkt-Annäherung der Bootstrap-Verteilung

Die Annäherung scheint ziemlich gut zu sein!

Wir könnten eine noch bessere Annäherung erhalten, indem wir die Sattelpunktannäherung und die Neuskalierung integrieren:

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

Die auf dieser Näherung basierende kumulative Verteilungsfunktion konnte nun durch numerische Integration gefunden werden, es ist jedoch auch möglich, eine direkte Sattelpunktnäherung dafür vorzunehmen. Aber das ist für einen anderen Beitrag, das ist lang genug.

Abschließend noch einige Kommentare aus der obigen Entwicklung. In wir eine Annäherung vorgenommen, bei der der dritte Term im Wesentlichen ignoriert wurde. Warum können wir das tun? Eine Beobachtung ist, dass für die normale Dichtefunktion der ausgelassene Term nichts beiträgt, so dass die Approximation genau ist. Da es sich bei der Sattelpunktnäherung um eine Verfeinerung des zentralen Grenzwertsatzes handelt, nähern wir uns dem Normalwert, sodass dies gut funktionieren sollte. Man kann sich auch konkrete Beispiele anschauen. Betrachtet man die Annäherung des Sattelpunkts an die Poisson-Verteilung und den ausgelassenen dritten Term, so wird dies in diesem Fall zu einer Trigammafunktion, die in der Tat eher flach ist, wenn das Argument nicht nahe bei Null liegt.(**)

Schließlich, warum der Name? Der Name stammt von einer alternativen Herleitung, die Techniken der komplexen Analyse verwendet. Später können wir das untersuchen, aber in einem anderen Beitrag!

kjetil b halvorsen
quelle
4
Was Sie bisher haben, ist großartig. Die Entwicklung dort ist sehr klar.
Glen_b
1
kjetil ich versucht , vier kleine Fehler zu beheben 1. „ Bei der Entwicklung wil ich folge “ 2 „ die für die approximatrion zur Arbeit “ 3. „ Was wir jetzt vermisst 4.“ „ implizite Differenzierung des sadlepoint “ aber in so tut es Es sieht so aus, als hätte ich eine Ihrer Gleichungen gebrochen. Ich habe keine Ahnung, wie, da ich nur diese Textelemente geändert habe (wie Sie aus dem Bearbeitungsverlauf ersehen können). Ich würde es zurücksetzen, aber da ich nicht erklären kann, wie das Beheben dieser Fehler zu einem Problem geführt hat, möchte ich keine weiteren Probleme verursachen. Entschuldigen Sie. (Es sah tatsächlich so aus, als wäre es kaputt
gegangen,
1
Möglicherweise gibt es einen mathJax-Fehler oder einen Fehler im Bearbeitungscode, der zu diesem Problem führt.
Glen_b
1
@Christoph Hanck: Um eine Annäherung an ein bestimmtes , lösen Sie die Sattelpunktgleichung , um zu finden . xt(§)t
kjetil b halvorsen
2
Vielleicht sollte darauf hingewiesen werden, dass die resultierende Sattelpunktnäherung außerhalb der konvexen Hülle der Daten undefiniert ist, wenn die empirische cgf verwendet wird. Siehe Feuerverger (1989) "On the Empirical Saddlepoint Approximation". Dies sollte auch im obigen Bootstrap-Beispiel der Fall sein.
Matteo Fasiolo
15

Hier werde ich auf die Antwort von kjetil eingehen und mich auf die Situationen konzentrieren, in denen die kumulatorenerzeugende Funktion (CGF) unbekannt ist, diese aber aus den Daten , wobei geschätzt werden kann . Der einfachste CGF-Schätzer ist wahrscheinlich der von Davison und Hinkley (1988) der in kjetils Bootstrap-Beispiel verwendet wird. Dieser Schätzer hat den Nachteil, dass die resultierende Sattelpunktgleichung nur dann gelöst werden kann, wenn , der Punkt, an dem die Sattelpunktdichte bewertet werden soll, in die konvexe Hülle von .x1,,xnxRd

K^(λ)=1ni=1neλTxi,
K^(λ)=y,
yx1,,xn

Wong (1992) und Fasiolo et al. (2016) haben dieses Problem gelöst, indem sie zwei alternative CGF-Schätzer vorgeschlagen haben, die so ausgelegt sind, dass die Sattelpunktgleichung für jedes gelöst werden kann . Die Lösung von Fasiolo et al. (2016), genannt Extended Empirical Saddlepoint Approximation ESA, ist im esaddle R-Paket implementiert , und ich gebe hier einige Beispiele.y

Betrachten Sie als einfaches univariates Beispiel die Verwendung von ESA, um eine Dichte von zu approximieren .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)

Das ist die Passform

Bildbeschreibung hier eingeben

Wenn wir den Teppich betrachten, ist es klar, dass wir die ESA-Dichte außerhalb des Bereichs der Daten bewertet haben. Ein schwierigeres Beispiel ist das folgende verzogene bivariate Gaußsche.

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

Bildbeschreibung hier eingeben

Die Passform ist ziemlich gut.

Matteo Fasiolo
quelle
9

Dank der großartigen Antwort von Kjetil versuche ich, selbst ein kleines Beispiel zu finden, über das ich gerne sprechen möchte, weil es einen relevanten Punkt aufzuwerfen scheint:

Betrachten Sie die -Verteilung. und seine Derivate werden gefunden hier und werden wiedergegeben in den Funktionen im Code unten. χ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)

Dies erzeugt

Bildbeschreibung hier eingeben

Dies ergibt offensichtlich eine Annäherung, die die qualitativen Merkmale der Dichte richtig macht, aber, wie in Kjetils Kommentar bestätigt, ist es keine richtige Dichte, da sie überall über der exakten Dichte liegt. Das erneute Skalieren der Approximation wie folgt ergibt den im Folgenden dargestellten, fast vernachlässigbaren Approximationsfehler.

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)

Bildbeschreibung hier eingeben

Christoph Hanck
quelle
1
Es ist ein Merkmal, dass die Sattelpunkt-Approximation nicht zu einem integriert werden muss, sondern oft nahe ist. Sie kann durch numerische Integration neu skaliert werden.
kjetil b halvorsen
Es könnte aufschlussreicher sein, den relativen Fehler aufzuzeichnen!
kjetil b halvorsen
approximationerror_unscaled/approximationerror_scaledEs stellt sich heraus, um 25.90798 zu schweben
Christoph Hanck