Greetings,
I've been playing with the umacs package for a few days and have worked out
an example of a simple linear regression using gibbs samplers (included
below). While extremely basic, I hope this might be helpful. I would love
to see more examples of MH sampling as well.
##############################
## Simple linear regression using UMacs package
##load libraries
library(Umacs);library(rv);library(MASS);library(mvtnorm)
##simulate data
n=500
b1 = 10
b2 = 6
beta = matrix(c(b1,b2),2,1)
sig = 4
x=cbind(1,rnorm(n,seq(1,20,length=n),sd=1))
y=matrix(rnorm(n,x%*%beta,sqrt(sig)),n,1)
#priors on beta and beta variance
b.prior = as.matrix(ncol=1,nrow=2,c(0,0))
v.prior = solve(diag(1000,2))
s1.prior = .1 # gamma prior
s2.prior = .1
#set up the initial values and define the dimension of each variable
b.init <- function () as.matrix(rnorm(2,0,1)) # b[1] is
intercept, b[2] is slope
v.init <- function () rnorm(1,0,1)^2 # error
#define the conditional probabilities for the gibbs sampler for the betas
and sigma
b.update <- function (){
#From Clark text "Models for Ecological Data" ISBN 0-691-12178-8 (page
192) and Clark workbook (pg 81)
# the If() statements allows this function to be used for a full
covariance matrix or just a single variance parameter
if(length(v) == 1) {
sx = crossprod(x) * 1/v
sy = crossprod(x,y) * 1/v
}
if(length(v) > 1) {
ss = t(x) %*% 1/v
sx = ss %*% x
sy = ss %*% y
}
bigv = solve(sx + v.prior)
smallv = sy + v.prior %*% b.prior
b = t(rmvnorm(1,bigv%*%smallv,bigv))
return(b)
}
v.update <- function () {
sx = crossprod((y - x%*%b))
u1 = s1.prior + 0.5*n
u2 = s2.prior + 0.5*sx
return(1/rgamma(1,u1,u2))
}
s <- Sampler(
.title = "Linear Regression",
x = x,
y = y,
b = Gibbs (b.update, b.init),
v = Gibbs (v.update, v.init),
Trace("b[1]"),
Trace("b[2]"),
Trace("v")
)
s(n.iter=2000)
#note original values are (approximately) returned
######################################################
Message: 48
Date: Wed, 31 Oct 2007 12:15:51 -0400
From: Ted Hart <[EMAIL PROTECTED]>
Subject: [R] Simple Umacs example help..
To: r-help@r-project.org
Message-ID: <[EMAIL PROTECTED]>
Content-Type: text/plain
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]]
--
Adam Wilson
http://hydrodictyon.eeb.uconn.edu/people/wilson/
Department of Ecology and Evolutionary Biology
BioPharm 223
University of Connecticut
Tel: 860.486.4157
[EMAIL PROTECTED]
[[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.
