21/03/12 14:23:51.96 TOSLyW/N.net
他でも使えるようにJAGS使って結果をグラフ表示できるように関数化しておいた。おもちゃ箱に入れておこうっと。
SideEffect <- function(p,n,k){
library(rjags)
library(BEST)
if(k/n > p){
lower=p
upper=1
}else{
lower=0
upper=p
}
'%&%' <- function(x,y) paste0(x,y)
modelstring='
model
{
theta ~ dunif(' %&% lower %&% ',' %&% upper %&% ')
k ~ dbin(theta,n)
}
'
writeLines(modelstring,'TEMPmodeljags.txt')
dataList=list(k=k,n=n)
jagsModel = jags.model( file="TEMPmodeljags.txt" ,data=dataList, quiet=FALSE)
update(jagsModel)
codaSamples=coda.samples(jagsModel,variable.names=c("theta"), n.iter=500000,thin=1)
js=as.matrix(codaSamples)
theta=js[,'theta']
plotPost(round(theta,7),xlab=quote(theta),col='lightseagreen')
invisible(theta)
}
SideEffect(0.05,10,3)