Chúng ta có thể sử dụng Base R để tìm 95% diện tích dưới đường cong không?

Question

Sử dụng Base R, tôi đã tự hỏi liệu tôi có thể xác định khu vực 95% dưới đường cong được biểu thị là posterior bên dưới không?

Cụ thể hơn, tôi muốn di chuyển từ mode (đường đứt nét màu xanh lá cây) về phía đuôi và sau đó dừng lại khi tôi đã che phủ 95% diện tích đường cong. Mong muốn là các giá trị trục x là giới hạn của khu vực 95% này như trong hình dưới đây?

 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 Nói cách khác, nó là rất hấp dẫn nếu một Interval như vậy là ngắn nhất khoảng 95% càng tốt.

Answer 1

Symmetric phân phối

Mặc dù dụ OP là không chính xác đối xứng, nó là đủ gần - và hữu ích để bắt đầu có kể từ khi giải pháp đơn giản hơn nhiều.

Bạn có thể sử dụng kết hợp integrate và optimize. Tôi đã viết điều này như là một chức năng tùy chỉnh, nhưng lưu ý rằng nếu bạn sử dụng điều này trong các tình huống khác, bạn có thể phải suy nghĩ lại các giới hạn để tìm kiếm số lượng.

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

Sử dụng nó như thế này,

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

cho

bất đối xứng phân phối

Trong trường hợp của một phân phối không đối xứng, chúng ta phải tìm kiếm hai điểm đó đáp ứng tiêu chí P (a < x < b) = Prob, trong đó Prob là xác suất mong muốn. Vì có nhiều khoảng thời gian vô hạn (a, b) thỏa mãn điều này, nên OP đề nghị tìm ra khoảng thời gian ngắn nhất.

Quan trọng trong giải pháp là định nghĩa domain, khu vực nơi chúng tôi muốn tìm kiếm (chúng tôi không thể sử dụng -Inf, Inf, vì vậy người dùng phải đặt giá trị này thành giá trị hợp lý).

# consider interval (a,b) on the x-axis 
# integrate our function, normalize to total area, to 
# get the total probability in the interval 
prob_ab <- function(fun, a, b, domain){ 
    totarea <- integrate(fun, domain[1], domain[2])[[1]] 
    integrate(fun, a, b)[[1]]/totarea 
} 

# now given a and the probability, invert to find b 
invert_prob_ab <- function(fun, a, prob, domain){ 

    O <- function(b, fun, a, prob){ 
    (prob_ab(fun, a, b, domain=domain) - prob)^2 
    } 

    b <- optimize(O, c(a, domain[2]), a = a, fun=fun, prob=prob)$minimum 

return(b) 
} 

# now find the shortest interval by varying a 
# Simplification: don't search past the mode, otherwise getting close 
# to the right-hand side of domain will give serious trouble! 
prob_int_shortest <- function(fun, prob, domain){ 

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

    # objective function to be minimized: the width of the interval 
    O <- function(a, fun, prob, domain){ 
    b <- invert_prob_ab(fun, a, prob, domain) 

    b - a 
    } 

    # shortest interval that meets criterium 
    abest <- optimize(O, c(0,mode), fun=fun, prob=prob, domain=domain)$minimum 

    # now return the interval 
    b <- invert_prob_ab(fun, abest, prob, domain) 

return(c(abest,b)) 
} 

Bây giờ, hãy sử dụng mã ở trên như sau. Tôi sử dụng một chức năng rất không đối xứng (chỉ giả sử mydist thực sự là một số pdf phức tạp, không phải là dgamma).

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) 

Trong ví dụ này, tôi đặt miền thành (0,10), vì rõ ràng khoảng thời gian phải ở trong một nơi nào đó. Lưu ý rằng việc sử dụng một giá trị rất lớn như (0, 1E05) không hoạt động, bởi vì integrate gặp sự cố với các chuỗi dài gần như số không. Một lần nữa, trong trường hợp của bạn, bạn sẽ phải điều chỉnh miền (trừ khi ai đó có ý tưởng tốt hơn!).

Answer 2

Dưới đây là một giải pháp làm cho việc sử dụng Trapezoidal rule.Bạn sẽ lưu ý rằng giải pháp được cung cấp bởi @Remko vượt trội hơn rất nhiều, tuy nhiên giải pháp này hy vọng thêm một số giá trị sư phạm khi nó chiếu sáng cách các vấn đề phức tạp có thể được giảm xuống thành hình học đơn giản, số học và các cấu trúc lập trình cơ bản như for loops.

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