# Example R code for rejection sampling, Metropolis sampling, and # Metropolis-Hastings sampling. # # M.V.Wickerhauser, 2011/03/02 /*************************************************** * Maxima code to find the max value of beta(2,2): * ***************************************************/ /* pdf without the constant */ f:x*(1-x); /* integrate to find the constant */ c:integrate(f,x,0,1); /* normalize the pdf */ f:(1/c)*f; /* differentiate to get f' */ diff(f,x); /* solve f'(x0)=0 to find max */ x0:solve(0=diff(f,x),x); /* evaluate f(x0) to get max value */ subst(x0,f); ############################################################## # From Seefeld and Linder, p.189 --- Rejection sampling, beta(2,2) rsamp<-function(n=5000, seed=NULL) { if( !is.null(seed) ) set.seed(seed); # in case a seed is provided x<-runif(n); # step 1 u<-runif(n); # step 2 fx<- dbeta(x,2,2); # step 3 acceptx<-x[fx>1.5*u]; # step 4 # scale by 1.5 = max(dbeta) hist(acceptx, breaks=20, prob=TRUE); points(x,fx); length(acceptx); # announce how many samples were accepted } ############################################################## # From Seefeld and Linder, p.203 --- Metropolis sampling, beta(4,3) nchain<-1000; y<-3; n<-5; theta<-vector(length=nchain); theta[1]<-0.04; alpha<-1+y; beta<-1+(n-y); # posterior shape parameters for(i in 2:nchain){ thetastar<-runif(1); # prior density is symmetric u<-runif(1); r<-thetastar**y*(1-thetastar)**(n-y)/(theta[i-1]**y*(1-theta[i-1])**(n-y)); ifelse( u