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 <ehart1@uvm.edu>
Subject: [R] Simple Umacs example help..
To: r-help@r-project.org
Message-ID: <A231C9C6-A711-4D0F-A2FB-1B634C802C42@uvm.edu>
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
Adam.Wilson@UConn.edu

	[[alternative HTML version deleted]]

