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);