BLOG

「役立てるだけの能力がなかった」

ベイズの定理

統計の基本定理の一つである「ベイズの定理」が発表されたのは1763年。いまではGoogle検索のアルゴリズムにも採用されるほど重宝されているが(検索履歴を事前分布として使っているものと思われる;詳細は不明)、実際に実用化が進んだのは早くとも2000年以降だと思われる。ベイズの定理が発表された当時は、膨大な計算量が必要とされることから、「役に立たない」ものとみなされていたらしい。ベイズの定理の後に、Fisherにより最尤法なる画期的な推定手法が発表されると、その簡便性と正確性の高さから広い支持を得た。長い長い間、ベイズの定理が表立って利用されることはほとんどなかった。

SPONSOR LINK

マルコフ連鎖モンテカルロ

しかし、ごく最近になって、マルコフ連鎖モンテカルロなる逐次的な計算手法が開発された。これは、ベイズの定理と至極相性がよい。というのも、ベイズの定理における計算量の問題が、この計算手法によってかなり緩和されるからだ(※)。もちろん、手計算で行える範疇ではないが、近年のPCの計算能力をもってすれば、十分計算できるほどのものになる。「役に立たない」と思われていたベイズの定理に、やっと人類が追いついた。計算量の問題さえ切抜けてしまえば、ベイズのその他もろもろの利点が浮き彫りになり、各種研究分野で実用化が進んでいる。Google検索におけるベイズ推定の応用もそのひとつになるのだろう。

「今」役に立つ

「役に立つ」という言葉は、今目に見えている問題を相手にとった言葉なのだろう。役に立つ研究というものが、今後起こりうる未知の問題に対し、どれだけの示唆を与えてくれるのだろうか、と思う(「役に立つ研究」は、今の問題には役に立つけど、それっきりなんだろうなぁ)。「質の高い、けど今は役に立たない基礎研究」をサポートする姿勢というは、将来起こりうる不測の事態に対する先行投資とも思うのだけれど。税金を使っている限り、それは過ぎたお願いだろうか。参照:「なんの役に立つんですか?の暴力性

※ ベイズの定理の計算には、パラメーターpが任意の値をとるすべての事象の数え上げが必要になる。このため、普通の計算手法では、PCをもってしても不可能なほどの計算が必要になる。

ナダルが全仏で10回優勝したすごさについて

ナダルがテニスの全仏オープンで10回目の優勝を達成しました。すごすぎて、このすごさを表現するのに適切な言葉が見つかりません。なので、確率で表してみたいと思います。

SPONSOR LINK

指数分布のFitting

全仏オープンといえば、テニスの4大大会のひとつに数えられる大会。出場する選手は全力でここにスケジュールを調整してくるし、基本的にトップ選手はみんなエントリーする。なので、世界で一番優勝の難しい大会の一つといえる。ナダルはその大会で10回目の優勝を達成した。これが全人未踏なのはもちろんなのだが、どれくらいの確率で起きうるのかが気になる。なんとか推定できないか。そこでウィキペディアをググると、なんと全仏オープンの歴代優勝者があるではないか。優勝経験者のみを対象に、ナダルの偉業がどれほど「偉業」なのか推定してみる。横軸にある個人の優勝回数をとり、縦軸にそれぞれの優勝回数を達成した人数をとる(ヒストグラム)。これに対して、指数分布をあてはめる;

f(x) = θexp(-θx)
(xは優勝回数)

どれくらいの確率でおきるのか

指数分布でFittingすると、θ = 0.55と推定され、以下のような図になる。縦軸は確率密度で、値が高いほど対応するx軸の事象が生じる可能性が高い。

よく聞く名前では、フェデラー、ジョコビッチは1回、ヴィランデルは3回、ボルグが6回優勝。ナダルの「10回」という優勝回数が生じる確率は、0.004、つまり0.4%だった。上の図でいうと、(ほとんど見えないけど)赤で塗りつぶされているところの面積と対応する。ちなみにこれは、「全仏優勝経験者」が10回優勝する確率。なので、1000年全仏オープンを続けたら、確率論的には4人10回優勝する人がでるかもね、っていうくらい(実際にはこれよりもずっと低いと思う;選手生命の長さなども考えれば)。地殻変動レベルの現象。これはすごい。ナダル、おめでとう。

負の二項分布モデルの精度検証:glmer.nb、glmmamb、JAGS

精度はいかほどですか

この前紹介したけれども、lme4パッケージの更新で、ランダム効果入りの負の二項分布モデルが扱えるようになっていた。しかし、他の関数(glmmadmbやJAGS)と比べてどの程度の精度なのか、気になるところ。そこで、簡単に推定の精度比較をしてみる。モデルは極めて単純な線形モデル。

Y ~ NB(pi, n)
pi = n/(mui + n)
log(mui) = a + bXi + εj
εj ~ Normal(0,σR2)

まず、上記を真のモデルとするシミュレーションデータを生成し、それをglmer.nbなどに推定させる、という作業を100回繰り返す。100回分の推定値の頻度分布と真の値を比べるということをした。パラメーターの真の値はそれぞれa=1.0, b=0.5, σR = 0.7とした。

結果

赤い曲線が100回分の推定値の頻度分布。縦の赤い実線は100回分の推定値の中央値、グレーの実線は真の値。左から切片(a)、回帰係数(b)、ランダム効果のばらつき(σR)。

glmer.nb

glmmadmb

JAGS

こうしてみると、推定精度の観点から一番難がありそうなのはglmmadmb。特にランダム効果のばらつきの推定値は、右側の裾が随分と重い。これはランダム効果が二つ以上になったら絶望的ではなかろうか。glmer.nbは比較的安定していたが、回帰係数の推定値のピークが平らになっているのがやや不安。中央値は真の値に近いけれども。JAGSは可もなく不可もなく、推定としては一番マシかな、という気がする。ランダム効果が増えたときが気になるけれども(たぶんJAGSの推定がかなりまマシなものになる)、面倒なのでやらない。でも、最尤推定の負の二項分布モデルで、ランダム効果の入れ過ぎはよくないということは間違いなさそう。確認してないけど。

シミュレーションデータ生成に使ったコード

glmer.nbの精度検証に使ったコード。そのほかは割愛。
apply系を使えばfor構文が要らなくなる気がするけど。。。

for(i in 1:100){
  #simulate 100 data points from 10 groups with one continuous fixed effect variable
  x<-runif(100,0,1)
  f<-gl(n = 10,k = 10)
  data<-data.frame(x=x,f=f)
  modmat<-model.matrix(~x,data)
  
  #the fixed effect coefficient
  p <- NULL
  p[1] <- 1; p[2] <- 0.5
  fixed<-c(p[1],p[2])
  
  #the random effect
  p[3] <- 0.7
  rnd<-rnorm(10,0,p[3])
  
  #the simulated response values
  size <- 3
  mu <- exp(modmat%*%fixed+rnd[f])
  prob <- size/(mu+size)
  data$y<-rnbinom(100,size=size,prob=prob)

  fit <- glmer.nb(y~x+(1|f),data)
  est[[1]][i] <- fit@beta[1]
  est[[2]][i] <- fit@beta[2]
  est[[3]][i] <- fit@theta
}

スーパー堤防

父の実家が陸前高田なので、節目で被災地に訪れている(特にボランティアというわけではないです)。そのときに聞いた興味深い話。

被災以降、津波の「減災」を目的としたスーパー堤防がガンガン作られており、それはもうすごい勢い。生態学の観点からそれが好ましくないことはもちろんそうなのだが、現地に生きる人の感情を考えると極めて難しい問題だなぁ、と感じていた。しかし、それは大きな勘違い(みんな堤防ほしいだろう、という先入観)で、現地からでる「堤防」に対する意見は住民の間でも大きく割れており、そしてそのパターンがどうにも奇妙だった。

奇妙だったというのは、海の近くに住んでいる住人(つまり被災リスクが高い)ほど堤防反対の意見をもつことが多かったこと。逆に、海から離れた場所に住む人(被災リスクは相対的に低いと思われる)ほど堤防建設を熱望しているという。

海際の人は年配の人が多く住んでいる可能性あるので、年齢層との偽相関は拭えないが、海際の住人の共通した意見として、海に住んでいるのだからそういったリスクは承知の上、海を見れないことのほうがよっぽどつらいし、仕事もしにくい、とのこと(もちろん、これが一般的な傾向というつもりはございません。あくまで聞いた範疇です)。

まったく定量的な調査ではないですが、現地の人に聞く限り、上述のような傾向があるそうです。だれかきちんと証明しているのかな、これ。

さまざまな情報量基準

生態学の界隈では、情報量基準なるものが大人気です。そのくせ、(筆者も含め)それぞれの情報量基準なるものがどんな意味をもつのか、よくわかっていない場合がほとんどです。わかっている範疇で整理。

AIC(赤池情報量基準)

情報量基準といったらAICと言わんばかりですが、なんだったら決定係数との違いすらよく分かっていないかったのが博士課程のころの僕です。AICが低いモデルでは、「同じ母集団から同じサンプル数の別のデータセットを取り出してきたときにも、それなりのパフォーマンス(当てはまり)を発揮する」ということになる。だが、これはサンプル数を無限大に増やしたとしても、必ずしも「真のモデル」には収束しないらしい。あくまで評価基準が予測精度だからだろうか。

AICc

サンプルサイズが小さいときはAICc!みたいな使われ方をしている。だが実は、正規分布を仮定した補正をしているので、正規分布以外のモデルに対して適用するのはご法度らしい。Burnhamさんの本を参照のこと。

BIC(ベイズ情報量基準)

AICと異なり、こちらは真のモデルを選ぶらしい(あくまでサンプルが無限大になったとき、ということだろうか)。あまり使ったことはないが、たしかにAICとは全く異なるモデルを選ぶ傾向がある。予測精度ではなく、真に重要な変数を選び出したい場合にはBICを使った方がよいのだろう。式的には、AICの罰則項にサンプルサイズが入っていないのに対し、BICには入っている。そのため、サンプルサイズが大きくなったときに、前者は対数尤度の項の影響が相対的に大きくなり、より複雑なモデルを選びがちになるようだ。

DIC(Deviance Information Criteria)

こちらはベイズモデルなどでよく使われる情報量基準。だが、広く使われている割に、その数学的根拠は随分と希薄らしい(論文でこういった表記をよく見かける)。DICの式の構造としては、モデルの当てはまりに関わる指標(Deviance)に対し、パラメーター数で罰則を与えている。AICでは2k(kはパラメーター数)という罰則項があるが、これには数学的根拠がある。DICの場合は、そういった裏付けがない、ということだろうか。時間のあるときに調べよう。
また、DIC利用の際の注意点として、パラメーターの正規性が担保された状態でしか役に立たないらしい。詳細はよくわからないが、ベイズのご利益がでる階層線形モデルや状態空間モデルでは、この「パラメーターの正規性」がほとんど担保されないらしい。パラメーターの正規性とはなんぞや。。。(おそらく、パラメーターの事後分布が正規分布で近似できる、ということか)。だとしたらDICを使える場面などそうそうないのでは。

wAIC(widely applicable information criteria)

最近考案されたらしいwAICなるものは、階層をもった複雑なモデルにおいても適用できるらしい。これまたAICに引き続き日本人が発見したので驚き。日本人って統計につおいのだろうか。開発者のワタナベさんのページはこちら(ワタナベさんが開発したからwAICかと思った。違うのかorz)。内容をまったく理解していないのですが、浅はかな理解では予測精度を問題にしているような気がするので、モデル選択の結果の解釈としてはAICと同様になるような気がする。階層モデルのようなものにも適用できる点が違いか(あくまでユーザー視点からすれば)。もう少し勉強する。

lme4パッケージの更新

しばらく見ないうちに、最新版のlme4パッケージでは劇的な進化を遂げていたようなので、それについてメモ。

glmer.nbの出現

これまで負の二項分布を仮定したモデルで、ランダム効果をいれるとしたら、glmmADMBパッケージか、もしくはベイズモデル化(WinBUGSなど)する必要があった(ADMBもベイズみたいなもん?中身はよくわからない)。しかし、最新版のlme4パッケージでは、glmer.nbという関数が表れたことにより、負の二項分布でもランダム効果を簡単に入れることができるようになった。仕様については基本的なlmerと変わらない模様。ただし、頻度主義的に扱うには尤度関数が複雑すぎる気がするので、いったいどこまで結果が信用できるかは不明。>3のランダム項をいれるのはきっとご法度。

mcmcsamp()の廃止

昔のlme4には、最尤法ではなく、MCMCでベイズっぽく推定するための関数が用意されていた。しかし、最新版のlme4パッケージでは除外された模様。Rのhelpには以下のように書いてある。
Due to difficulty in constructing a version of mcmcsamp that was reliable even in cases where the estimated random effect variances were near zero (e.g. https://stat.ethz.ch/pipermail/r-sig-mixed-models/2009q4/003115.html), mcmcsamp has been withdrawn (or more precisely, not updated to work with lme4 versions >=1.0.0)

どうやらランダム効果のばらつきがゼロ付近のときに、あまりいい推定をしてくれないということで、やめにしたということらしい。個人的には、ランダム効果が多いとき(>3)の最尤推定は不安なので、mcmcsamp()の取り下げは手痛い、、、JAGSとかはコード書くのが面倒だからなぁ、、、

lmerにおけるREMLのデフォルト化

REMLとはなんぞや、と調べてみると、REstricted Maximum Likelihoodの略らしい。最新版のlmerでは、MLではなくREMLがデフォルトになっている。REMLとMLの違いは、(筆者の浅はかな理解によると)、前者は固定効果による自由度減少を考慮したうえでランダム効果の分散を推定し、後者は考慮しないらしい。つまり、後者の場合、ランダム効果のばらつきの推定は偏りが生じているということになる(不偏推定量でない?)。ただ、REMLではモデル選択ができないため(なぜかはよくわからない)、MuMInパッケージなどでモデル選択する場合はREML=Fとする必要がある。正規分布以外ではMLしか使えないので心配の必要なし。

glmerとlmerの(よくわからない)分離

正規分布の場合はlmer、それ以外はglmerということでよさそう。

QGIS: 属性テーブルの一括編集

QGISで属性テーブルの値を一括編集したいときの処理

1.変更対象の地物を選択する
2.属性テーブルからフィールド演算をひらく
3.既存カラムを更新を選択
4.変更したいカラムを選択
5.演算スペースに一括入力したい文字や数値を入力。
※文字の場合、シングルコーテーション(’ ‘)で囲む
※選択フィーチャーだけに適用のチェックボックスにチェック

ATMの悪夢

住み慣れた国を離れると、これまでの常識が通じないことも多々あり、いろいろと苦労しております。2017年の4月からアメリカに移動し、所属がUniversity of Minnesotaとなったのですが、さっそく新生活の立ち上げでひやり(というか実際損した)とした経験です。

日本ではお金を振り込んだり、預け入れをする際に、ATMの機械がその場で紙幣を数えてくれるのが普通だと思います。しかし、アメリカではそうではありませんでした。Depositのボタンを押すと、なにやら封筒がいるか、と聞いてくるのです。こちらの常識で考えると、当然封筒はいらないので、いらないと押し、そのまま紙幣を挿入口へ。ただウイーンといいながら吸い込まれていく紙幣。まったく数える気配がありません。よくわからないので、そのままその場を後にしたのですが、そもそも紙幣を数えていなかったこと、そして挿入金額が自己申告(途中で入力させられる)だったことに違和感を覚えるべきでした。不安になり、しばらく経ってからググってみると、、、げっ、、、封筒に金額と名前を記入して入れるだと、、、?オンラインアカウントの口座残高をみると、入金されていなことになっている(ちなみに500ドル)。なんとアメリカでは、ATMの預け入れは夜間金庫みたいなものになっており、行員があとでATMからお金を運び出し、手作業で数えているとのこと。なんじゃそりゃ。最近では日本のようなATMが大多数のようですが、運悪く私の使ったATMは旧式のもので、従来の方式が採用されていた模様。

そんな馬鹿な、と思い、支店にいって事情を話すと、行員の方が丁寧に対応してくれました。しかし時すでに遅し、紙幣をそのまま挿入したためか、500ドル入金したはずが480ドルになってしまいました。とほほ。全部なくなるよりはマシか。。。

発見率モデル(除去法)

除去法による発見率モデル

発見率を推定する方法として、繰り返し調査の他に、除去法による推定方法があります。除去法では、調査区間から個体を一時的に除去していき(例えば、バケツに入れておくなど)、個体数の減少パターンから発見率を推定します。個体は閉鎖区間からいずれいなくなるのですが(全数調査:この場合推定はいらない)、それを実行するのは現実的ではありません。そのため、2回~3回の調査を行い、その限られた情報から真の個体数を推定することになります。
 除去法で得られたデータは、次のような過程を経て観察されていると考えることができます(簡単のため、3回除去を行った場合を想定)。まず、3回調査したときにとれた総個体数nは、真の個体数Nに採捕効率π(0–1の値をとる:真の個体数に対する採捕個体数の割合)をかけたものなので、

n = πN

と表すことができます。次に、1回目の調査で採捕された個体数n1は、真の個体数Nに対し、「調査一回あたりの採捕効率p」を掛け合わせたものになります。

n1 = Np = nπ−1 p

2回目の調査で採捕される個体数n2は、真の個体数Nから1回目の調査で採捕された個体数を差し引いたもの(N− n1)に対し、採捕効率pを掛け合わせたものとして表します。

n2 = (N− n1)p = (nπ−1 − nπ−1 p)p = nπ−1p(1−p)

同様に考えて、3回目の調査で採捕される個体数n3は以下の式で表せます。

n3 = (N− n1− n2)p = nπ−1p(1−p) (1−p)

つまり、k回目の採捕個体数を一般化した形で表すと、

nk = nπ−1p(1−p)k−1

となります。

このように得られるデータは、多項分布(発見プロセス)とポアソン分布(真の個体数Nの空間的ばらつき)の混合モデルとして表現することができます。

[ni1, ni2, ni3|Ni, π] ~ Multinomial(Ni, π)

Ni ~ Pois(λ)

BUGSによる実装

実際にスクリプトを書く場合、BUGSの多項分布は限られた観察数n(総採捕数)に対する各要素の比率(nk/n)を推定することしかできないので、真の個体数Nと総採捕個体数nの間の関係を記述する際に、二項分布を仮定することになります。

model{
###Priors###------------------------------------------------
p0 ~ dbeta(1,1)
log(lambda) <- alpha
alpha ~ dnorm(0,ninfo)
  
###Likelihood###--------------------------------------------
  for(i in 1:Nsite){
    #p[i]: probablity of being captured with a single pass
    p[i] <- p0
   
    mu[i,1] <- p[i] # proportion of fish removed at 1st pass
    mu[i,2] <- p[i]*(1-p[i])
    mu[i,3] <- p[i]*(1-p[i])*(1-p[i])
   
    #total catchability
    pi[i] <- sum(mu[i,])
    
    for(j in 1:J){# J = 3
      muc[i,j]	<- mu[i,j]/pi[i]
    }
    y[i,1:J] ~ dmulti(muc[i,1:J], ncap[i])
    ncap[i] ~ dbin(pi[i], N[i]) #ncap: total number of fish captured
  }
  
  for(i in 1:Nsite){ N[i] ~ dpois(lambda) }
}

発見率について、地点間でランダムにばらつく場合を想定する場合、以下のように書き換えます。

model{
###Priors###------------------------------------------------
logit.p0 <- logit(p0)
p0 ~ dbeta(1,1)
log(lambda) <- alpha
alpha ~ dnorm(0,ninfo)
tau.eps <- pow(sigma.eps,-2)
sigma.eps ~ dunif(0,100)

###Likelihood###--------------------------------------------
  for(i in 1:Nsite){
    #p[i]: probablity of being captured with a single pass
    logit(p[i]) <- logit.p0 + eps[i]
    eps[i] ~ dnorm(0,tau.eps)
   
    mu[i,1] <- p[i] # proportion of fish removed at 1st pass
    mu[i,2] <- p[i]*(1-p[i])
    mu[i,3] <- p[i]*(1-p[i])*(1-p[i])
   
    #total catchability
    pi[i] <- sum(mu[i,])
    
    for(j in 1:J){# J = 3
      muc[i,j]	<- mu[i,j]/pi[i]
    }
    y[i,1:J] ~ dmulti(muc[i,1:J], ncap[i])
    ncap[i] ~ dbin(pi[i], N[i]) #ncap: total number of fish captured
  }
  
  for(i in 1:Nsite){ N[i] ~ dpois(lambda) }
}

そのほか、説明変数に関する情報があれば、それらの発見率に対する効果なども見ることが可能です。

発見率モデル(繰り返し法)

ロバストデザインによる調査

生物の分布や個体数を調べ、対象生物のある地点における存在確率を推定する際、「いない」ことを証明するのはとても難しい問題です。なぜなら、生物調査で「本当はいる個体」を100%発見することは不可能だからです。この問題をそのままにすると存在確率を過小評価し、パラメーター推定に偏りをもたらすことが知られています。

この問題に対応するため、最近ではロバストデザインによる調査が一般的になってきています。ロバストデザインとは、同じ調査地点で複数回調査する調査デザインのことで、これにより発見率に起因する「在・不在」や「個体数」データのばらつきを統計的に補正することができます。基本的な原理としては、繰り返し調査の間は「真の在・不在」の状態は変化しないという仮定を置くことで、「同じ地点における調査回ごとの在・不在データのばらつきは発見率に由来する」というモデリングを行います。

BUGSによるモデリング(応答変数が在・不在)

調査回tにおける調査地点iでの在・不在をYti、真の在不在をzi(調査回tの間ではzは変化しないと仮定するので、tは省く)、発見率をpとしたとき、在・不在データが得られる確率は下記のように整理されます。

真の状態は在、データも在:zi*p
真の状態は在、データは不在:zi*(1-p)
真の状態は不在、データも不在:zi*p or zi*(1-p)

これらをモデル化するためには、まずは真の状態(zi)を記述し、次に真の状態ziに依存して観察値Ytiが得られるプロセスをモデル化します。

model{
#Likelihood
for(i in 1:Nsite){
 Y[i] ~ dbin(mu[i],T.survey)  #Y[i]はT.survey回調査したうち、在となった回数
 mu[i] <- p*z[i] #データが在となる確率は真の状態と発見率の積
 z[i] ~ dbern(psi) #真の存在確率psi
}
# Prior
p ~ dunif(0,1)
psi ~ dunif(0,1)

}

BUGSモデルでは、Ytiは発見回数として表されます(例えば、3回調査して2回見つかった場合はY=2)。前半の二項モデルの部分(Y[i] ~ dbin(mu[i],T.survey))では、発見率に応じて3回中何回が「在」となるそうかをモデル化しています。

真の在・不在の状態に影響する要因を調べたい場合は、下記のようにスクリプトを書き替えることで容易にパラメーター推定を行えます。

model{
#Likelihood
for(i in 1:Nsite){
 Y[i] ~ dbin(mu[i],T.survey)
 mu[i] <- p*z[i]

 z[i] ~ dbern(psi[i])
 logit(psi[i]) <- b[1] + b[2]*x[i]
}

# Prior
p ~ dunif(0,1)
for(i in 1:2){ b[i] ~ dnorm(0, 0.001) }

}

なお、実際にモデルを回す際には、zの初期値の与え方に注意しなければなりません。zがありえない値(例えばY=3なのにz=0となるなど)になるとBUGSはエラーコードを返してきます。そのため、zの初期値はz = rep(1, N.sample)などとしておくほうが無難です。

BUGSによるモデリング(応答変数が個体数)

応答変数が個体数などの離散値の場合でも、少しスクリプトを書き替えるだけで対応することが可能です。調査回tにおける調査地点iでの発見個体数をYti、真の在不在をNiとすると:

model{
#Likelihood
for(i in 1:Nsite){
 for(t in 1:Nsurvey){
  Y[t,i] ~ dbin(p, N[i])
 }
 N[i] ~ dpois(lambda) #真の個体数の平均値lambda
}

# Prior
p ~ dunif(0,1)
log(lambda) <- log.lambda
log.lambda ~ dnorm(0, 0.001)
}

なお、ここでの発見率pは、発見された個体数を真の個体数で除したもの(見つかった個体の割合)として解釈することができます。上と同様に、lambdaを線形予測子と関連付けることで、真の個体数に影響する要因を調べることができます。

発見率のばらつき

上記のモデルでは、場所によらず発見率は一定としていましたが、様々な理由で発見率が地点によってばらつくことが考えられます。そのような場合には、発見率をなんらかの確率分布に従うランダム変数とすることで、そのばらつきを考慮することが可能です。上記の在・不在モデルにおいて、発見率がランダムにばらつくと仮定したモデルは次のように書くことができます。

model{
#Likelihood
for(i in 1:Nsite){
 Y[i] ~ dbin(mu[i],T.survey)
 mu[i] <- p[i]*z[i]
 z[i] ~ dbern(psi)

 logit(p[i]) <- logit.mean.p + eps[i]
 eps[i] ~ dnorm(0, tau.p)
}
# Prior
logit(mean.p) <- logit.mean.p
mean.p ~ dunif(0,1)
psi ~ dunif(0,1)
tau.p ~ dgamma(0.001,0.001)

}

上のモデルではロジットスケールで正規分布を仮定した誤差を考えていますが、ベータ分布によるモデリングなども可能です。