関数解読discint

Bayesian Computation with R: Second Edition (Use R!)

離散確率分布に対して最高確率区間を計算する。

　最高確率区間：こちらから与える確率を含む最短区間
　　確率分布が単峰でない場合は，2つ以上の区間に分かれる。
　　確率分布が左右対称でない場合は，区間の位置が分布中心から左右のどちらかにずれる。

使用法

discint(dist, prob)

引数
dist:1列目に確率変数の値，2列目に対応する確率を含む行列
prob:最高確率区間を求めるためにこちらから与える確率の値

値
prob:下記のsetに含まれる確率変数のどれかが実現する正確な確率（>=prob）
set:最高確率区間に含まれる確率変数の値の集合

discint = function (dist, prob)
#関数discintの中身
{

    x = dist[, 1]
    p = dist[, 2]
     # distの中身をx,pに分けて保存。
    n = length(x)
     # xの長さをnに保存。
     
    sp = sort(p, index.return = TRUE)
     # pを昇順で並べ替えて，対応する要素の位置も出力し，
     #  spにリストとして保存。

    ps = sp$x
     # 昇順に並べ替えたpをpsに保存。後で使わず。
     
    i = sp$ix[seq(n, 1, -1)]
     # seq(n, 1, -1)
     #  nから1まで1ずつ減少する数列
     # sp$ix[seq(n, 1, -1)]
     #  pの昇順に対応する添字を降順に変換しiに保存。
     # 結局，iにはpの要素を降順で並べ替えた場合の
     #  pの要素の位置が保存されている。


    ps = p[i]
     # pの要素を降順で並べ替えてpsに保存。
     # ps = sp$x[seq(n,1,-1)]でも良さそう。

    xs = x[i]
     # xをpの要素を降順に対応するように並べ替えてxsに保存。

    cp = cumsum(ps)
     # psの累積和を計算。
     #  cp[i]=sum(ps[1:i])

    ii = 1:n
     # 別の添字1〜nの作成。

    j = ii[cp >= prob]
     # cp>=probを満たすcpの要素について，iiの位置でTRUEと表す。

    j = j[1]
     # cp>=probを最初に満たしたcpの要素の位置をjに保存。
     #  cp[j]>=prob かつ cp[j-1]<prob が成立。

    eprob = cp[j]
     # cp[j]をeprobに保存。

    set = sort(xs[1:j])
     # xs[1:j]を昇順に並べ替え，setに保存。

    v = list(prob = eprob, set = set)
     # eprobとsetをリスト形式にしてvに保存。
    
    return(v)
     # vを出力。
}

例

x=0:10
probs=dbinom(x,size=10,prob=.3)

# 図示
plot(x,probs,type="h")

dist=cbind(x,probs)
pcontent=.8

(hpi = discint(dist,pcontent))
$prob
[1] 0.8214841

$set
[1] 1 2 3 4