Hello all...
I am just starting to teach myself Bayesian methods, and am
interested in learning how to use UMacs. I've read the
documentation, but the single example is a bit over my head at the
level I am at right now. I was wondering if anyone has any simple
examples they'd like to share. I've successfully done a couple of
simple gibbs examples, but have had a hard time modifying some of the
home written metropolis hastings samplers I've made to work with
Umacs. Does anyone have any pointers or simple 2 parameter
examples? Thanks.
Here is one of my simple MH samplers using a simple linear regression
with a Cauchy error term.
x <- c(1.808,1.799,1.179,0.574,3.727,0.946,3.719,1.566,3.596,3.253)
y <- c(1.816,1.281,-1.382,0.573,3.793,0.935,1.775,1.474,3.679,3.889)
fn = function(x,a=0,b=1){
a+b*x
}
sample.ab <-function(x,y,a,b,s,da,db){
bstar = runif(1,b-db,b+db)
astar = runif(1,a-da,a+da)
logalpha = sum(dcauchy(y,location=fn(x,astar,bstar),scale=s,log=T) -
dcauchy(y,location=fn(x,a,b),scale=s,log=T))
logu = log(runif(1,0,1))
acc = (logu < logalpha)
b = acc*bstar + (1-acc)*b
a = acc*astar + (1-acc)*a
list(b=b,a=a,acc=acc)
}
samples = function(x,y,a,b,s,ds){
sstar = runif(1,s-ds,s+ds)
while(sstar <= 0){
sstar = runif(1,s-ds,s+ds)
}
logalpha = sum( dcauchy(y,location=fn(x,a,b),scale=sstar,log=T) -
dcauchy(y,location=fn(x,a,b),scale=s,log=T)) - log(sstar/s)
logu = log(runif(1,0,1))
acc = (logu < logalpha)
s = acc*sstar + (1-acc)*s
list(s=s,accs=acc)
}
sample.abs<-function(n=10000,x,y,a=0,b=1,s=2,da=.2,db =.2,ds=1)
{
accab <- 0
accs <- 0
A = B = S = rep(NaN,n)
for(i in 1:n){
z = sampleab(x,y,a,b,s,da,db)
q <- samples(x,y,a,b,s,ds)
A[i] = a = z$a
B[i] = b = z$b
S[i]=s=q$s
accab = accab + z$acc
accs <- accs +q$accs
}
invisible(list(a=A, b=B, s=S, accab=accab/n,accs=accs/n))
}
Cheers,
Ted
Dept. of Biology,
University of Vermont
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
