## The file for two examples in Part 1 # Do cauchyprior() # Do nmusigma.gibbs() f <- function(x) { u <- (x-0.0675)^2 exp(-10 * u) /(1+x^2) } f1 <- function(x) { u <- (x-0.0675)^2 x * exp(-10 * u) /(1+x^2) } canal <- function() { bot <- integrate(f, lower=-Inf, upper=Inf) top <- integrate(f1, lower=-Inf, upper=Inf) ex <- top$integral/bot$integral ex } cauchyprior <- function(N=10000, tstart=0, n=20, tmu=0, tsd=1){ set.seed(44) dat <- rnorm(n, mean=tmu, sd = tsd) # we simulate n data points # if you have your own data just pass it to this routine x <- rep(NA, N) x[1] <- tstart ybar <- mean(dat) cat("Sample mean=", round(ybar, 4), "\n") y <- rcauchy(N) for (i in 2:N) { top <- dnorm(y[i], mean=ybar, sd=1/sqrt(n)) # Numerator in the ratio bot <- dnorm(x[i-1], mean=ybar, sd=1/sqrt(n)) # Denominator in the ratio alpha <- min(1, top/bot) # Acceptance probability if (runif(n=1) < alpha) # Test whether to accept x[i] <- y[i] else x[i] <- x[i-1] } u <- post.proc(x) tmean <- canal() z <- list(tmean=tmean, gmean=u[1], nse=u[2], lcor=u[3]) z } nmusigma.gibbs <- function(N=2000, tstart=0, s2start=1.0, n=20, tmu=5){ set.seed(44) dat <- rnorm(n, mean=tmu) # we simulate n data points # if you have your own data just pass it to this routine theta <- rep(NA, N) s2 <- rep(NA, N) xbar <- mean(dat) sy2 <- sum((dat-xbar)^2) z <- matrix(NA, 2, 4) dimnames(z) <- list(c("mu", "sigma2"), c("tmean", "gmean", "nse", "lag-1.cor") ) z[1,1] <- xbar #var(theta) <- sy2/(n(n-3)) z[2, 1] <- sy2/(n-3) #var (sigma2) <- 2 * z[2, 1]^2/(n-5) s2[1] <- s2start theta[1] <- tstart cat("Sample mean=", round(xbar, 4), "\n") for (i in 2:N) { theta[i] <- rnorm(n=1, mean=xbar, sd=sqrt(s2[i-1]/n)) u <- sum( (dat - theta[i])^2 ) /2 s2[i] <- 1.0/rgamma(n=1, shape=n/2, rate=u) } u <- post.proc(theta) v <- post.proc(s2) z[1, 2:4] <- u z[2, 2:4] <- v z <- round(z, 4) #list(theta=theta, s2=s2, dat=dat) z } post.proc <- function(x, nburn = length(x)/4) { # x is the MCMC sample # nburn is the number of initial iterations to discard N <- length(x) z <- x[nburn:N] print(summary(z)) #Obtain summary par(mfrow = c(1, 2)) tsplot(z) # Check the time series plot v <- acf(z, lag.max = 50, plot = T) u <- calcbatch2(z, 20) par(mfrow = c(1, 1)) u[-4] } calcbatch2 _ function(a, batchsiz) { # # calcbatch2 -- standard errors by batch means method # # Author: Kate Cowles # l <- length(a) k <- batchsiz m <- floor(l/k) # number of batches meanvect <- numeric() for(i in 0:(m - 1)) { meanvect[i + 1] <- mean(a[(k * i + 1):(k * (i + 1))]) } myvar <- numeric() myvar[1] <- mean(meanvect) myvar[2] <- sqrt(var(meanvect)/m) myvar[3] <- acf(meanvect, 1, "corr", F)$acf[2] myvar[4] <- k names(myvar) <- c("Mean", "SEM", "Lag 1 autocorr between batches", "Batch size") myvar }