Bagaimana cara pendekatan saddlepoint bekerja?

38

Bagaimana cara pendekatan saddlepoint bekerja? Masalah apa yang cocok untuk itu?
(Jangan ragu untuk menggunakan contoh atau contoh tertentu dengan cara ilustrasi)

Apakah ada kekurangan, kesulitan, hal-hal yang harus diperhatikan, atau perangkap untuk yang tidak waspada?

Glen_b -Reinstate Monica
sumber

Jawaban:

49

Perkiraan saddlepoint ke fungsi densitas probabilitas (berfungsi juga untuk fungsi massa, tapi saya hanya akan berbicara di sini dalam hal kepadatan) adalah perkiraan yang bekerja dengan sangat baik, yang dapat dilihat sebagai penyempurnaan pada teorema limit pusat. Jadi, itu hanya akan bekerja di pengaturan di mana ada teorema batas pusat, tetapi perlu asumsi yang lebih kuat.

Kita mulai dengan asumsi bahwa fungsi pembangkit momen ada dan dua kali dapat dibedakan. Ini menyiratkan secara khusus bahwa semua momen ada. Biarkan menjadi variabel acak dengan fungsi pembangkit momen (mgf) dan cgf (fungsi pembangkit kumulans) (di mana menunjukkan logaritma natural). Dalam pengembangan saya akan mengikuti erat Ronald W Butler: "Saddlepoint Approximations with Applications" (CUP). Kami akan mengembangkan pendekatan saddlepoint menggunakan pendekatan Laplace ke integral tertentu. Menulis X

M(t)=EetX
K(t)=logM(t)log
eK(t)=etxf(x)dx=exp(tx+logf(x))dx=exp(h(t,x))dx
mana . Sekarang kita akan memperluas Taylor di mempertimbangkan sebagai konstanta. Ini memberi mana ' menunjukkan diferensiasi sehubungan dengan x . Perhatikan bahwa h '(t, x) = - t- \ frac {\ partial} {\ partial x} \ log f (x) \\ h' '(t, x) = - \ frac {\ partial ^ 2} {\ partial x ^ 2} \ log f (x)> 0 (ketidaksetaraan terakhir dengan asumsi karena diperlukan agar aproksimasi berfungsi). Biarkan 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
xtmenjadi solusi untuk h(t,xt)=0 . Kami akan menganggap bahwa ini memberikan minimum untuk h(t,x) sebagai fungsi dari x . Dengan menggunakan ekspansi ini di bagian integral dan lupa tentang bagian , beri
eK(t)exp(h(t,xt)12h(t,xt)(xxt)2)dx=eh(t,xt)e12h(t,xt)(xxt)2dx
yang merupakan integral Gaussian, memberikan
eK(t)eh(t,xt)2πh(t,xt).
Ini memberikan (versi pertama) dari perkiraan saddlepoint sebagai
(*)f(xt)h(t,xt)2πexp(K(t)txt)
Perhatikan bahwa aproksimasi memiliki bentuk keluarga eksponensial.

Sekarang kita perlu melakukan beberapa pekerjaan untuk mendapatkan ini dalam bentuk yang lebih bermanfaat.

Dari kita mendapatkan Membedakan ini sehubungan dengan memberikan (dengan asumsi kami), jadi hubungan antara dan adalah monoton, jadi didefinisikan dengan baik. Kita membutuhkan perkiraan untuk . Untuk itu, kita dapatkan dengan menyelesaikan darih(t,xt)=0

t=xtlogf(xt).
xt
txt=2xt2logf(xt)>0
txtxtxtlogf(xt) log f ( x t ) = K ( t ) - t x t - 1(*)
(**)logf(xt)=K(t)txt12log2π2xt2logf(xt).
xtxtlogf(xt) Dengan asumsi istilah terakhir di atas hanya bergantung pada , sehingga turunannya terhadap kira-kira nol (kami akan kembali untuk mengomentari ini), kami mendapatkan Sampai dengan perkiraan ini maka kita memiliki sehingga dan harus dihubungkan melalui persamaan yang disebut persamaan saddlepoint. xtxt
logf(xt)xt(K(t)xt)txtt
0t+logf(xt)xt=(K(t)xt)txt
txt
(§)K(t)xt=0,

Apa yang kita lewatkan sekarang dalam menentukan adalah dan dapat kita temukan dengan diferensiasi implisit dari persamaan saddlepoint : Hasilnya adalah bahwa (hingga perkiraan kami) Menyatukan semuanya, kami memiliki perkiraan akhir saddlepoint untuk kepadatan sebagai h ( t , x t ) = - 2 log f ( x t )(*)

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).
Sekarang, untuk menggunakan ini secara praktis, untuk memperkirakan kerapatan pada titik tertentu , kita memecahkan persamaan saddlepoint untuk untuk menemukan .xtxtt

The saddlepoint pendekatan sering dinyatakan sebagai pendekatan dengan kepadatan rata-rata berdasarkan pengamatan iid . Fungsi penghasil kumulans dari mean adalah hanya , jadi pendekatan saddlepoint untuk mean menjadi nX1,X2,,XnnK(t)

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

Mari kita lihat contoh pertama. Apa yang kita dapatkan jika kita mencoba mendekati kerapatan normal standar mgf adalah jadi sehingga persamaan saddlepoint adalah dan pendekatan saddlepoint memberikan jadi dalam hal ini aproksimasi tepat.

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

Mari kita lihat aplikasi yang sangat berbeda: Bootstrap di domain transformasi, kita dapat melakukan bootstrap secara analitik menggunakan pendekatan saddlepoint ke distribusi bootstrap rata-rata!

Asumsikan kita memiliki iid didistribusikan dari beberapa kepadatan (dalam contoh simulasi kita akan menggunakan distribusi satuan eksponensial). Dari sampel kami menghitung fungsi penghasil momen empiris dan kemudian cgf empirik . Kita membutuhkan mgf empiris untuk rata-rata yaitu dan cgf empiris untuk mean yang kita gunakan untuk membangun pendekatan saddlepoint. Dalam beberapa kode R berikut (R versi 3.2.3): X1,X2,,Xnf

M^(t)=1ni=1netxi
K (t)=log M (t)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

(Saya sudah mencoba menulis ini sebagai kode umum yang dapat dimodifikasi dengan mudah untuk cgf lain, tetapi kodenya masih belum terlalu kuat ...)

Kemudian kami menggunakan ini untuk sampel sepuluh pengamatan independen dari unit distribusi eksponensial. Kami melakukan bootstrap nonparametrik yang biasa "dengan tangan", plot histogram bootstrap yang dihasilkan untuk rata-rata, dan overplot perkiraan pendekatan 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)

Memberikan plot yang dihasilkan:

pendekatan saddlepoint distribusi bootstrap

Perkiraannya agak bagus!

Kita bisa mendapatkan perkiraan yang lebih baik dengan mengintegrasikan pendekatan saddlepoint dan men-rescaling:

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

Sekarang fungsi distribusi kumulatif berdasarkan perkiraan ini dapat ditemukan dengan integrasi numerik, tetapi juga dimungkinkan untuk membuat perkiraan saddlepoint langsung untuk itu. Tapi itu untuk posting lain, ini sudah cukup lama.

Akhirnya, beberapa komentar tidak disertakan dalam pengembangan di atas. Dalam kami melakukan perkiraan yang pada dasarnya mengabaikan istilah ketiga. Kenapa kita bisa melakukan itu? Satu pengamatan adalah bahwa untuk fungsi kepadatan normal, istilah yang ditinggalkan tidak berkontribusi apa-apa, sehingga perkiraannya tepat. Jadi, karena saddlepoint-aproksimasi adalah penyempurnaan pada teorema limit pusat, jadi kita agak mendekati normal, jadi ini seharusnya bekerja dengan baik. Kita juga dapat melihat contoh-contoh spesifik. Melihat pendekatan saddlepoint ke distribusi Poisson, melihat istilah ketiga yang ditinggalkan itu, dalam hal ini yang menjadi fungsi trigamma, yang memang agak datar ketika argumennya tidak mendekati nol.(**)

Akhirnya, mengapa namanya? Nama ini berasal dari derivasi alternatif, menggunakan teknik analisis kompleks. Nanti kita bisa melihat itu, tetapi di pos lain!

kjetil b halvorsen
sumber
4
Apa yang Anda miliki sejauh ini hebat. Perkembangan di sana sangat jelas.
Glen_b -Reinstate Monica
1
kjetil Saya mencoba untuk memperbaiki empat kesalahan ketik kecil 1. " Dalam pengembangan saya akan mengikuti " 2. " diperlukan untuk perkiraan untuk bekerja " 3. " Apa yang kita lewatkan sekarang " 4. " diferensiasi implisit dari sadlepoint " tetapi dalam melakukannya sepertinya saya melanggar salah satu persamaan Anda - Saya tidak tahu caranya, karena saya tidak mengubah apa pun selain item teks tersebut (seperti yang dapat Anda lihat dari riwayat edit). Saya akan mengembalikannya tetapi karena saya tidak bisa menjelaskan bagaimana memperbaiki kesalahan-kesalahan itu menyebabkan masalah, saya tidak ingin menyebabkan masalah lebih lanjut. Permintaan maaf saya.
1
Mungkin saja ada bug mathJax atau bug dalam kode edit yang mengarah ke masalah ini.
Glen_b -Reinstate Monica
1
@Christoph Hanck: Untuk mendapatkan perkiraan di xif , Anda memecahkan persamaan untuk menemukan . xt(§)t
kjetil b halvorsen
2
Mungkin perlu menunjukkan bahwa, ketika cgf empiris digunakan, perkiraan saddlepoint yang dihasilkan tidak terdefinisi di luar lambung cembung data. Lihat Feuerverger (1989) "Pada Pendekatan Saddlepoint Empiris". Ini juga harus terjadi pada contoh bootstrap di atas.
Matteo Fasiolo
15

Di sini saya memperluas jawaban kjetil, dan saya fokus pada situasi-situasi di mana Cumulant Generating Function (CGF) tidak diketahui, tetapi dapat diperkirakan dari data , di mana . Estimator CGF paling sederhana mungkin adalah dari Davison and Hinkley (1988) yang digunakan dalam contoh bootstrap kjetil. Estimator ini memiliki kelemahan bahwa persamaan saddlepoint yang dihasilkan dapat diselesaikan hanya jika , titik di mana kami ingin mengevaluasi kepadatan saddlepoint, termasuk dalam cembung cembung . x R d K ( λ ) = 1x1,,xnxRd K '(λ)=y,yx1,...,xn

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

Wong (1992) dan Fasiolo et al. (2016) mengatasi masalah ini dengan mengusulkan dua alternatif penduga CGF, yang dirancang sedemikian rupa sehingga persamaan saddlepoint dapat diselesaikan untuk setiap . Solusi dari Fasiolo et al. (2016), disebut ESA Pendekatan Empiris Saddlepoint diperpanjang, diimplementasikan dalam paket esaddle R dan di sini saya memberikan beberapa contoh.y

Sebagai contoh univariat yang sederhana, pertimbangkan untuk menggunakan ESA untuk memperkirakan kepadatan .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)

Ini cocok

masukkan deskripsi gambar di sini

Melihat permadani itu jelas bahwa kami mengevaluasi kepadatan ESA di luar rentang data. Contoh yang lebih menantang adalah Gaussian bivariat berikut yang melengkung.

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

masukkan deskripsi gambar di sini

Cocok cukup bagus.

Matteo Fasiolo
sumber
9

Terima kasih atas jawaban Kjetil yang luar biasa, saya mencoba memberikan sedikit contoh sendiri, yang ingin saya diskusikan karena sepertinya memunculkan poin yang relevan:

Pertimbangkan . dan turunannya dapat ditemukan di sini dan direproduksi dalam fungsi dalam kode di bawah ini. K ( t )χ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)

Ini menghasilkan

masukkan deskripsi gambar di sini

Ini jelas menghasilkan perkiraan yang mendapatkan fitur kualitatif dari kepadatan yang tepat, tetapi, seperti yang dikonfirmasi dalam komentar Kjetil, bukan kepadatan yang tepat, karena di atas kepadatan yang tepat di mana-mana. Menyesuaikan ulang perkiraan sebagai berikut memberikan kesalahan perkiraan yang hampir diabaikan yang digambarkan di bawah ini.

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)

masukkan deskripsi gambar di sini

Christoph Hanck
sumber
1
Ini adalah fitur, pendekatan saddlepoint tidak perlu diintegrasikan ke satu, tetapi sering dekat. Itu dapat ditingkatkan dengan integrasi numerik.
kjetil b halvorsen
Itu bisa lebih mengungkapkan untuk merencanakan kesalahan relatif!
kjetil b halvorsen
approximationerror_unscaled/approximationerror_scaledternyata berkisar sekitar 25,90798
Christoph Hanck