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
XM(t)=EetX
K(t)=logM(t)logeK(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)=−tx−logf(x)h(t,x)xth(t,x)=h(t,x0)+h′(t,x0)(x−x0)+12h′′(t,x0)(x−x0)2+⋯
′xh′(t,x)=−t−∂∂xlogf(x)h′′(t,x)=−∂2∂x2logf(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)(x−xt)2)dx=e−h(t,xt)∫∞−∞e−12h′′(t,xt)(x−xt)2dx
yang merupakan integral Gaussian, memberikan
eK(t)≈e−h(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)=0t=−∂∂xtlogf(xt).
xt∂t∂xt=−∂2∂x2tlogf(xt)>0
txtxt∂∂xtlogf(xt) log f ( x t ) = K ( t ) - t x t - 1(*)
logf(xt)=K(t)−txt−12log2π−∂2∂x2tlogf(xt).(**)
xtxt∂logf(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)∂t∂xt−t
0≈t+∂logf(xt)∂xt=(K′(t)−xt)∂t∂xt
txtK′(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)∂x2t=−∂∂xt(∂logf(xt)∂xt)=−∂∂xt(−t)=(∂xt∂t)−1
K′(t)=xt∂xt∂t=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π−−√e−12x2
M(t)=exp(12t2)K(t)=12t2K′(t)=tK′′(t)=1
t=xtf(xt)≈e12t2−txt12π⋅1−−−−−√=12π−−√e−12x2t
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,…,XnfM^(t)=1n∑i=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:
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!
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,…,xn x∈Rd
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)
Ini cocok
Melihat permadani itu jelas bahwa kami mengevaluasi kepadatan ESA di luar rentang data. Contoh yang lebih menantang adalah Gaussian bivariat berikut yang melengkung.
Cocok cukup bagus.
sumber
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)
Ini menghasilkan
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.
sumber
approximationerror_unscaled/approximationerror_scaled
ternyata berkisar sekitar 25,90798