Jak działa przybliżanie saddlepoint?

Odpowiedzi:

49

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)=EetX
i cgf (funkcja generowania skumulowanego) K(t)=logM(t) (gdzie logoznacza 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
eK(t)=etxf(x)dx=exp(tx+logf(x))dx=exp(h(t,x))dx
gdzie h(t,x)=txlogf(x) . Teraz Taylor rozszerzyh(t,x) wx traktująct jako stałą. Daje to
h(t,x)=h(t,x0)+h(t,x0)(xx0)+12h(t,x0)(xx0)2+
gdzie Oznacza różnicowanie w odniesieniu dox . Zauważ, że
h(t,x)=txlogf(x)h(t,x)=2x2logf(x)>0
(ostatnia nierówność z założenia, ponieważ jest to potrzebne do przybliżenia, aby zadziałało). Niechxtbędzie rozwiązaniem dlah(t,xt)=0. Zakładamy, że daje to minimum dla h(t,x)w funkcjix. Używając tego rozszerzenia w całce i zapominając oczęści, daje
eK(t)exp(h(t,xt)12h(t,xt)(xxt)2)dx=eh(t,xt)e12h(t,xt)(xxt)2dx
która jest całką Gaussa, dającą
eK(t)eh(t,xt)2πh(t,xt).
Daje to (pierwsza wersja) przybliżenie punktu siodłowego jako
(*)f(xt)h(t,xt)2πexp(K(t)txt)
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,xt)=0 otrzymujemy

t=xtlogf(xt).
Zróżnicowanie tego względemxtdaje
txt=2xt2logf(xt)>0
(o przyjętych założeń), a więc stosunek międzytixtznaczy monotoniczne takxtjest dobrze zdefiniowana. Potrzebujemy przybliżenia doxtlogf(xt). W tym celu otrzymujemy rozwiązanie z log f ( x t ) = K ( t ) - t x t - 1(*)
(**)logf(xt)=K(t)txt12log2π2xt2logf(xt).
Zakładając, że powyższy ostatni termin słabo zależy tylko odxt, więc jego pochodna względemxtwynosi w przybliżeniu zero (wrócimy do komentowania tego), otrzymujemy
logf(xt)xt(K(t)xt)txtt
Do tego przybliżenia mamy wtedy
0t+logf(xt)xt=(K(t)xt)txt
, abytixtbyły powiązane przez równanie
(§)K(t)xt=0,
który nazywa się równaniem saddlepoint.

(*)

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

nX1,X2,,XnnK(t)

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

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

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!

X1,X2,,Xnf

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

saddlepoint przybliżenie dystrybucji bootstrap

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.

(**)

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
4
To, co do tej pory masz, jest świetne. Rozwój tam jest bardzo wyraźny.
Glen_b
1
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
1
Możliwe, że w kodzie edycji występuje błąd mathJax lub błąd, który prowadzi do tego problemu.
Glen_b
1
xt(§)t
2
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
15

x1,,xnxRd

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

y

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

wprowadź opis zdjęcia tutaj

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)

wprowadź opis zdjęcia tutaj

Dopasowanie jest całkiem dobre.

Matteo Fasiolo
źródło
9

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ę:

χ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

wprowadź opis zdjęcia tutaj

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)

wprowadź opis zdjęcia tutaj

Christoph Hanck
źródło
1
Jest to cecha, przybliżenie punktu saddlepoint nie musi być zintegrowane z jednym, ale często jest bliskie. Można go przeskalować za pomocą integracji numerycznej.
kjetil b halvorsen
Bardziej odkrywcze może być wykreślenie błędu względnego!
kjetil b halvorsen
approximationerror_unscaled/approximationerror_scaledOkazuje się, że unosi się około 25.90798
Christoph Hanck