ベースRを使ってカーブの下の面積の95％を見つけることはできますか？

Question

ベースRを使用して、下にposteriorと表示された曲線の下の95％の面積を判断できるかどうか疑問に思っていましたか？

具体的には、mode（緑色の破線）からテールに向かって移動し、カーブ領域の95％をカバーしたときに停止します。下の図に示すように、この95％の範囲の限界であるx軸の値が望ましいですか？

 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ような間隔が可能インターバル最短95％である場合、それは非常に望ましいです。解決策は非常に簡単であることから、そこを開始するのに便利 -

Answer 1

OPの例では、正確に対称ではなかったにもかかわらず対称分布が

、それは十分に近いです。

integrateとoptimizeの組み合わせを使用できます。これをカスタム関数として書きましたが、これを他の状況で使用すると、分位点を検索するための境界を再考する必要があるかもしれないことに注意してください。

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

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

が

に非対称分布の場合には非対称分布

を与え、このようにそれを使用して、我々はその2点を検索する必要がP（a < x < b）= Prob、ここで、Probはある望ましい確率である。これを満たす間隔は無限に多いので（a、b）、OPは最短のものを見つけるよう提案した。

domainは、私たちが検索したい領域です（-Inf, Infは使用できませんので、ユーザーはこれを妥当な値に設定する必要があります）。

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

ここでこのようなコードを使用してください。私は非常に非対称な関数を使用しています（mydistは実際にはdgammaではなく複雑なpdfです）。

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) 

この例では、ドメインを（0,10）に設定していますが、明確に間隔がどこかにある必要があるためです。 integrateはゼロに近い長いシーケンスに問題があるため、（0,1E05）のような非常に大きな値を使用すると機能しないことに注意してください。あなたの状況に合わせて、ドメインを調整する必要があります（誰かが良いアイデアを持っていない限り！）。ここで

Answer 2

Trapezoidal ruleを利用したソリューションです。@Remkoによって提供されるソリューションははるかに優れていることに気づくでしょうが、このソリューションは、複雑な問題を単純なジオメトリ、算術、および基本的なプログラミング構造（例えば、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