関数解読pdisc

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

離散事前分布を用いた場合の比率に関する事後分布の計算

使用法

pdisc(p, prior, data)

引数
p:（推測する）比率の候補となる値で構成されるベクトル
prior:候補となる各比率の事前確率で構成されるベクトル
data:成功と失敗の数で構成されるベクトル
　確率pで成功が起こるベルヌーイ分布でモデル化

値
候補となる各比率の事前確率で構成されるベクトル

#関数pdiscの中身
pdisc = function (p, prior, data)
{
    s = data[1]
    f = data[2]
    p1 = p + 0.5 * (p == 0) - 0.5 * (p == 1)
     #ベクトルpの要素に0があれば，その要素に0.5を加える。
     #ベクトルpの要素に1があれば，その要素から0.5を引く。
     # 上の操作は何を意味する？
     #log(0)は計算できないからか。
     #ベクトルpの要素に0も1もなければ，p1=pでこれが通常。

    like = s * log(p1) + f * log(1 - p1)
     #対数尤度(A)   （定数項を除く）
       #確率p1で成功がs回発生
       #確率1-p1で失敗がf回発生

    like = like * (p > 0) * (p < 1) -
    999 * ((p == 0) * (s > 0) + (p == 1) * (f > 0))
     #通常は第1項目のみでOK。
     # p>0かつp<1の場合は対数尤度(A)をそのまま代入。

     # p==0かつs>0 （成功の確率は0なのに成功した結果あり）
     # p==1かつf>0 （失敗の確率は0なのに失敗した結果あり）
     #  この場合，log(0)=-Infの代わりに-999を代入。

    like = exp(like - max(like))

     #対数尤度を尤度に戻す。
     # max(like)を引く意味は？
     # 後でsum(post)=1となるよう基準化するので，
     # -max(like)はなくても同じ。
     # exp(like - max(like))=exp(like)*exp(-max(like))
     #  exp()の値を発散させないよう念のため付けてる？

    product = like * prior
     #事後分布のカーネル

    post = product/sum(product)
     #事後確率
     #sum(post)=1がとなるようにproductを基準化


    return(post)#事後確率postを出力
}

例

#比率（成功確率）と対応する事前確率
p=c(.0,.25,.3,1)
prior=c(.25,.25,.25,.25)

#成功と失敗の数
data=c(5,10)

#比率（成功確率）の事後確率
pdisc(p,prior,data)