Eğrinin altındaki alanın% 95'ini bulmak için Base R kullanabilir miyiz?

Question

Taban R'yi kullanarak, aşağıdaki posterior olarak belirtilen eğrinin altındaki% 95 alanını belirleyip belirlemediğimi merak ediyordum?

Daha ayrıntılı olarak, mode'dan (yeşil kesik çizgi) kuyruklara doğru gitmek istiyorum ve eğri alanının% 95'ini kapladığımda durmak istiyorum. Aşağıdaki resimde gösterildiği gibi bu% 95 alanın sınırları olan x ekseni değerleri istenir mi? Böyle bir aralık kısa% 95 aralığı mümkün Diğer bir deyişle

 prior = function(x) dbeta(x, 15.566, 7.051) 
likelihood = function(x) dbinom(55, 100, x) 
posterior = function(x) prior(x)*likelihood(x) 

mode = optimize(posterior, interval = c(0, 1), maximum = TRUE, tol = 1e-12)[[1]] 

curve(posterior, n = 1e4) 

P.S., oldukça tercih edilir. çözüm daha basit olduğu için ve orada başlatmak için faydalıdır -

Answer 1

OP'ın örnek tam olarak simetrik değildi bile Simetrik dağılım

, yeterince yakın.

integrate ve optimize'un bir kombinasyonunu kullanabilirsiniz. Bunu özel bir işlev olarak yazdım, ancak bunu başka durumlarda kullanırsanız, quantile'ı aramak için sınırları yeniden düşünmeniz gerekebilir.

# For a distribution with a single peak, find the symmetric! 
# interval that contains probs probability. Search over 'range'. 
f_quan <- function(fun, probs, range=c(0,1)){ 

    mode <- optimize(fun, interval = range, maximum = TRUE, tol = 1e-12)[[1]] 

    total_area <- integrate(fun, range[1], range[2])[[1]] 

    O <- function(d){ 
    parea <- integrate(fun, mode-d, mode+d)[[1]]/total_area 
    (probs - parea)^2 
    } 
    # Bounds for searching may need some adjustment depending on the problem! 
    o <- optimize(O, c(0,range[2]/2 - 1E-02))[[1]] 

return(c(mode-o, mode+o)) 
} 

böyle kullanın

,

f <- f_quan(posterior, 0.95) 
curve(posterior, n = 1e4) 
abline(v=f, col="blue", lwd=2, lty=3) 

asimetrik dağılım durumunda

Asimetrik dağılımını

verir, biz iki puan aramak zorunda P (a) < x < b) = Prob, Prob istenen bir olasılıktır. Sonunda çok fazla aralık olduğu için (a, b), OP en kısa olanı bulmayı önerdi.

Çözümde önemli olan, domain tanımının bulunduğu bölgedir, (-Inf, Inf'u kullanamıyoruz, böylece kullanıcı bunu makul değerlere ayarlamak zorundayız). Yukarıdaki kodu yukarıdaki gibi kullanın. Çok asimetrik bir fonksiyon kullanıyorum (sadece mydist'in aslında karmaşık bir pdf olduğunu, dgamma'yı değil).

mydist <- function(x)dgamma(x, shape=2) 
curve(mydist(x), from=0, to=10) 
abline(v=prob_int_shortest(mydist, 0.9, c(0,10)), lty=3, col="blue", lwd=2) 

Bu örnekte, aralığın açıkça bir yerde olması gerektiğinden, etki alanını (0,10) olarak ayarlıyorum. (0, 1E05) gibi çok büyük bir değerin kullanılmadığını unutmayın, çünkü integrate, sıfıra yakın uzun dizilerde sorun yaşar. Yine, durumunuz için, alanı düzeltmeniz gerekecek (birisi daha iyi bir fikri olmadığı sürece!).

Answer 2

İşte Trapezoidal rule bir çözüm yapım kullanımıdır.@Remko tarafından sağlanan çözümün çok daha üstün olduğunu, ancak bu çözümün, karmaşık problemlerin basit geometri, aritmetik ve for loops gibi temel programlama yapılarına ne kadar azaltılabileceğini aydınlattığı için bazı pedagojik değerlere umuyoruz.

findXVals <- function(lim, p) { 
    ## (1/p) is the precision 

    ## area of a trapezoid 
    trapez <- function(h1, h2, w) {(h1 + h2) * w/2} 

    yVals <- posterior((1:(p - 1))/p) 
    m <- which.max(yVals) 
    nZ <- which(yVals > 1/p) 

    b <- m + 1 
    e <- m - 1 
    a <- f <- m 

    area <- 0 
    myRng <- 1:(length(nZ)-1) 
    totArea <- sum(trapez(yVals[nZ[myRng]], yVals[nZ[myRng+1]], 1/p)) 
    targetArea <- totArea * lim 

    while (area < targetArea) { 
     area <- area + trapez(yVals[a], yVals[b], 1/p) + trapez(yVals[e], yVals[f], 1/p) 
     a <- b 
     b <- b + 1 
     f <- e 
     e <- e - 1 
    } 

    c((a - 1)/p, (f + 1)/p) 
} 

findXVals(.95, 10^5) 
[1] 0.66375 0.48975