STAT 705, Fall 2010 MC vs MCMC: Some Examples in R ================================= 0. Standard Normal and Double Exponential: General Information ================================================================ Double Exponential (DE) ------------------------ f(x)=0.5*exp(|-x|), -Inf < x < Inf F(x)=0.5*exp(x), x < 0 =0.5+0.5*(1-exp(-x)) x >= 0 #Note: E(X)=0, Var(X)=2 x <- seq(-5,5,0.01) fDE <- function(x){0.5*exp(- abs(x))} #Double Exponential pdf & cdf FDE <- function(x){ifelse((x <0), y <- 0.5*exp(x), y <- 0.5+0.5*(1-exp(-x)))} Standard Normal (SN) --------------------- fSN <- function(x){(1/sqrt(2*pi))*exp(-0.5*x^2)} #Standard Normal pdf ##Plots par(mfrow=c(2,2)); par(oma=c(0,0,4,0)) Plot of SN vs DE pdf's ----------------------- x <- seq(-4,4,0.01) plot(x,fDE(x),type="l", ylab="",lty=1) lines(x,fSN(x),type="l", lty=2) legend(-4.7,0.45,c("DE","SN"), lty=1:2, bty="n") title("SN & DE pdf's") Plot of SN vs DE cdf's ----------------------- plot(x,FDE(x),type="l", ylab="",lty=1) lines(x,pnorm(x),type="l", lty=2) legend(-4.7,0.85,c("DE","SN"), lty=1:2, bty="n") title("SN & DE cdf's") Plot Ratio fSN(x)/fDE(x) ------------------------- Ratio <- function(x) {fSN(x)/fDE(x)} plot(x,Ratio(x), type="l", main="fSN/fDE") #USA (in S-Plus): usa(states=T, coast=T, add=F, xlim=c(65,135), ylim=c(24,50), fifty=F) > mtext("Standard Normal (SN) vs Double Exponential (DE)", cex=1.5, side=3,outer=T, line=1) ###FIG 1 -------- I. MC (Rejection Method) for Generating N(0,1) Data From Centered DE ===================================================================== @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ fSN Note: c >= sup ---- = 1.315489 fDE Note: First generate Z = F^{-1}(U) ~ DE Note: Use ifelse((U <= f(Z)/(c*h(Z))), Z, NA) then eliminate NA's GenerateSNfromDE <- function(n,c){ #Get at most n N(0,1). c >=1.315489 #Generate DE Z from F^{-1}(U) U <- runif(n) Z <- ifelse((0 < U) & (U <0.5), log(2*U), -log(1-2*(U-0.5))) #Rejection method f <- function(x){(1/sqrt(2*pi))*exp(-0.5*x^2)} h <- function(x){0.5*exp(- abs(x))} U <- runif(n) Y <- ifelse((U <= f(Z)/(c*h(Z))), Z, NA) ###c >= 1.315489 Y <- Y[!is.na(Y)] #Get rid of NA's. Y ~ N(0,1) Y} @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ > Y <- GenerateSNfromDE(1000,5) unix.time(Y <- GenerateSNfromDE(1000,5)) user system elapsed 0.004 0.000 0.008 > length(Y) [1] 205 #For c=5, only 205 good N(0,1)'s out of 1000 iterations. hist(Y, prob=T, main="c=5: Generated Y. 205/1000") lines(x,fSN(x), type="l") #Check by qqplot: qqnorm(Y, main="c=5. Y Normal ?") abline(0,1) ##Now use smaller c: c=1.35 to get many more N(0,1)'s Y <- GenerateSNfromDE(1000,1.35) length(Y) [1] 750 hist(Y, prob=T, main="c=1.35: Generated Y. 750/1000") lines(x,fSN(x), type="l") #Check by qqplot: qqnorm(Y, main="c=1.35. Y Normal ?") abline(0,1) mtext("Rejection Method: N(0,1) From Double Exponential", cex=1.5, side=3,outer=T, line=1) ### FIG 2 --------- II. MCMC for Generating N(0,1) Data From Shifted DE ==================================================== Note: Pi is the stationary distribution from which we need data Note: Cannot use q()!!! Instead denote the proposal dist by q=qq Note: Proposal q=qq dist is shifted DE Note: qq(x,y)=qq(y,x); this gives Metropolis @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Pi <- function(x){(1/sqrt(2*pi))*exp(-0.5*x^2)} #Stationary dist. qq <- function(x,y){0.5*exp(-abs(x-y))} #Proppsal: shifted DE a <- function(x,y,Pi,qq){pmin(1,(Pi(y)*qq(x,y))/(Pi(x)*qq(y,x)))} SNfromDE.MCMC <- function(n){ X <- rep(0,n) for(t in 1:(n-1)){ U <- runif(1) #Generate Shifted DE Y from F^{-1}(U) Y <- ifelse((0 < U) & (U <0.5), log(2*U*exp(X[t])), -log(2*(1-U)*exp(-X[t]))) X[t+1] <- ifelse((runif(1) <= a(X[t],Y,Pi,qq)), Y, X[t])} X} @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ X <- SNfromDE.MCMC(5000) unix.time(X <- SNfromDE.MCMC(5000)) user system elapsed 1.552 0.005 1.570 hist(X, probability=T,main=""); title("Hist of X[1:5000]") z <- seq(-3.2,3.2, length=500) lines(z,Pi(z),type="l") #TS Plots ts.plot(X[2001:2300], main="X[2001:2300]") abline(2,0,lty=2) abline(-2,0,lty=2) qqnorm(X, main="X[1:5000]") abline(0,1) acf(X,main="Correlated Data!") ##Autocorrelation of X(t) mtext("MCMC: N(0,1) From Shifted Double Exponential. n=5000",line=1, side = 3, outer = T, cex = 1.5) ### FIG 3 --------- #Take every 10th X_t Note: To reduce the correlation we sample every 10th X_t XX <- rep(0,500) for(k in 1:500){XX[k] <- X[10*k]} hist(XX, probability=T,main=""); title("XX(k)=X[10*k]") z <- seq(-3.2,3.2, length=500) lines(z,Pi(z),type="l") #TS Plots ts.plot(XX[101:300], main="XX[101:300]") abline(2,0,lty=2) abline(-2,0,lty=2) qqnorm(XX, main="Every 10th X(t)") abline(0,1) acf(XX) ##Autocorrelation of XX(t) mtext("MCMC: N(0,1) From Shifted DE. Every 10th X(t)",line=1, side = 3, outer = T, cex = 1.5) ### FIG 4 --------- III. MCMC for Generating N(0,1) Data From N(X,S^2) ===================================================== par(mfrow=c(3,2)); par(oma=c(0,0,4,0)) @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ Note: Similar to II, but here the proposal q=qq is N(X,S^2) Note: Again Metropolis: qq(x,y)=qq(y,x) ss <- 0.5 #ss=variance Pi <- function(x){(1/sqrt(2*pi))*exp(-0.5*x^2)} #Stationary dist. qq <- function(x,y){(1/(sqrt(2*pi*ss)))*exp(-0.5*(x-y)^2/ss)} #Proppsal a <- function(x,y,Pi,qq){pmin(1,(Pi(y)*qq(x,y))/(Pi(x)*qq(y,x)))} SNfromN.MCMC <- function(n){ X <- rep(0,n) for(t in 1:(n-1)){ Y <- rnorm(1, X[t], sqrt(ss)) X[t+1] <- ifelse((runif(1) <= a(X[t],Y,Pi,qq)), Y, X[t])} X} @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ X <- SNfromN.MCMC(5000) > unix.time(X <- SNfromN.MCMC(5000)) user system elapsed 1.334 0.000 1.364 #TS Plots ts.plot(X[1001:1300], main="X[2001:2300], ss=0.5") abline(2,0,lty=2) abline(-2,0,lty=2) hist(X, probability=T, main=""); title("Hist of X[1:5000]") z <- seq(-3.2,3.2, length=400) lines(z,Pi(z),type="l") qqnorm(X, main="X[1:5000], ss=0.5") abline(0,1) acf(X, main="Correlated X") #Take every 10th X_t XX <- rep(0,500) for(k in 1:500){XX[k] <- X[10*k]} qqnorm(XX, main="Every 10th X(t)") abline(0,1) acf(X) acf(XX) mtext("MCMC: N(0,1) From N(X,S^2=0.5)",line=1, side = 3, outer = T, cex = 1.5) ### FIG 5 --------- %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% IV. Application: MCMC for E[h(W)], W Truncated Normal ====================================================== a. True -------- X ~ N(0,1), W= X| 0 < X < 1 meanW <- (1-exp(-0.5))/((sqrt(2*pi))*(pnorm(1)-0.5)) #True E[W] > meanW [1] 0.4598622 #True E[W] b. Check E[W] from data ------------------------ x <- rnorm(2000) w <- x[(x > 0) & (x < 1)] #Truncation mean(w) > length(w) [1] 671 ##Since there are less w's than x's > mean(w) [1] 0.4580269 #OK. Close to E[W]=0.4598622 c. Comparison of Histograms & pdf's ------------------------------------ Note: Comparison of the pdf's of N(0,1) and N(0,1) truncated to (0,1). par(oma = c(0, 0, 5, 0)) par(mfrow = c(1, 2)) hist(x, probability=T,main =""); title("X~N(0,1)") z <- seq(-3.2,3.2, length=500) f <- function(z) exp(-z^2/2)/sqrt(2*pi) lines(z,f(z),type="l") hist(w, probability=T,main =""); title("W ~ X | 0 0) & (w < 1))*exp(-w^2/2)/(sqrt(2*pi)*(pnorm(1)-0.5))} a <- function(w,y,Pi){pmin(1,Pi(y)/Pi(w))} TruncN.MCMC <- function(n,s){ ##Need n and s=SD W <- rep(0.1,n) for(t in 1:(n-1)){ Y <- rnorm(1, W[t],s) #Sample from Q(.|W[t])~N(W[t],s) W[t+1] <- ifelse((runif(1) <= a(W[t],Y,Pi)), Y, W[t])} W} @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ##Several cases for different s=SD par(mfrow = c(2,2)); par(oma=c(0,0,4,0)) ##s=0.5: W <- TruncN.MCMC(3000,0.5) > mean(W[1001:3000]) > [1] 0.4621351 #True E[W]=0.4598622 ts.plot(W[1001:1200], type="l") abline(0.45986,0,lty=2) title("SD=0.5. AVE[1001:3000]=0.4621351") unix.time(W <- TruncN.MCMC(3000,0.5)) user system elapsed 0.789 0.001 0.790 ##s=1: W <- TruncN.MCMC(3000,1) > mean(W[1001:3000]) [1] 1] 0.4606069 #True E[W]=0.4598622 ts.plot(W[1001:1200], type="l") abline(0.45986,0,lty=2) title("s=1. AVE[1001:3000]=0.4606069") ##s=5: W <- TruncN.MCMC(3000,5) > mean(W[1001:3000]) [1] 0.4577798 #True E[W]=0.4598622 ts.plot(W[1001:1200], type="l") abline(0.45986,0,lty=2) title("SD=5. AVE[1001:3000]=0.4577798") ##s=8: W <- TruncN.MCMC(3000,8) > mean(W[1001:3000]) [1] 0.3658715 #True E[W]=0.4598622 ts.plot(W[1001:1200], type="l") abline(0.45986,0,lty=2) title("SD=8. AVE[1001:3000]=0.3658715") mtext("200 iterations from Metropolis algorithm for stationary dist truncated N(0,1). E[W]=0.45986",line=1, side = 3, outer = T, cex = 1.5) ### FIG 7 ---------- #Hist & Acf: hist(W[1:3000],probability=T, main="Histogram of W[1:3000]") zt <- seq(0,1, length=500) f1 <- function(zt) exp(-zt^2/2)/(sqrt(2*pi)*(pnorm(1)-0.5)) lines(zt,f1(zt),type="l") hist(W[1001:3000],probability=T, main="Histogram of W[1001:3000]") lines(zt,f1(zt),type="l") hist(W[2001:3000],probability=T, main="Histogram of W[1001:3000]") lines(zt,f1(zt),type="l") acf(W) #Autocor. of W[t], t=1,...,3000 mtext("W(t) from Metropolis algorithm: Q~N(W(t),1), Stationary dist truncated N(0,1). E[W]=0.45986",line=2, side = 3, outer = T, cex = 1.5) ### FIG 8 ---------- e. Direct method ------------------ U <- runif(3000) W <- qnorm(0.5 + U*(pnorm(1)-0.5)) mean(W) > mean(W) [1] 0.459125 #Precise and fast!!! unix.time(W <- qnorm(0.5 + U*(pnorm(1)-0.5))) user system elapsed 0.001 0.000 0.001 ##Very fast. ##Hist & Acf: par(mfrow=c(2,1)); par(oma=c(0,0,7,0)) hist(W,probability=T, main="Histogram of W[1:3000]") zt <- seq(0,1, length=400) f1 <- function(zt) exp(-zt^2/2)/(sqrt(2*pi)*(pnorm(1)-0.5)) lines(zt,f1(zt),type="l") acf(W) mtext("W = X | 0 < X < 1, X~N(0,1), E[W]=0.45986, from a direct method: W <- qnorm(0.5 + U*(pnorm(1)-0.5))" ,line=2, side= 3, outer = T, cex = 1.5) ### FIG 9 ----------