Example R programs and commands Gibbs Sampling, finite discrete and bivariate normal # All lines preceded by the "#" character are my comments. # All other left-justified lines are my input. ######### Finite Discrete joint, marginal and conditional pdfs ## See the textbook, p.193, for an explanation of this method f<- matrix( c(0.3,0.2,0.1,0.2,0.1,0.1),2,3,byrow=TRUE); # joint pdf fA <- colSums(f); # marginal pdf for A fB <- rowSums(f); # marginal pdf for B fAgB <- diag(1/fB) %*% f; # conditional pdf for A given B fBgA <- f %*% diag(1/fA); # conditional pdf for B fiven A A<- c(1,2,3); # Prob. space for a B<- c(1,2); # Prob space for b a<-1; b<-1; # initialize (a,b) samples<-matrix(0, 2,3); # accumulate samples here for(i in 1:10000) { a<- sample(A,1, prob=fAgB[b,]); # replace a b<- sample(B,1, prob=fBgA[,a]); # replace b samples[b,a] <- samples[b,a] + 1; # count this (a,b) } ######### Bivariate normal joint and conditional pdfs ## See the textbook, p.195, for an explanation of this method gibbsBVN <- function(x,y, n=1000, rho) { m<-matrix(0, ncol=2,nrow=n); m[1,]<-c(x,y); # store initial values for(i in 2:n) { x<- rnorm(1, rho*y, sqrt(1-rho*rho)); y<- rnorm(1, rho*x, sqrt(1-rho*rho)); m[i,]<-c(x,y); # store new values } m } sim<-gibbsBVN(0,0,1000,rho=-0.9); plot(sim); # Alternative: use mvrnorm() require(MASS); rho<- -0.9; sig<-matrix(c(1,rho,rho,1),2,2); sim2<- mvrnorm(1000, mu=c(0,0), Sigma=sig) plot(sim2);