\documentclass[11pt]{article}

\usepackage{indentfirst}
\usepackage[utf8]{inputenc}

\RequirePackage{amsmath}
\DeclareMathOperator{\var}{var}

\begin{document}

\title{Stat 8931 Spin Glass Homework Solution, Part II}
\author{Charles J. Geyer}
\maketitle

\section{Setup}

Read betas.
<<read-write-beta>>=
foo <- try(scan("betas.txt"))
if (inherits(foo, "try-error")) {
    write(c(br, bd), file = "betas.txt")
    foo <- scan("betas.txt")
}
n <- sqrt(length(foo) / 2)
# n
br <- matrix(foo[1:n^2], n, n)
bd <- matrix(foo[n^2 + 1:n^2], n, n)
# br
# bd
@

Set tau.
<<tau-one>>=
tau <- 0.2
@

Setup.
<<setup>>=
i <- matrix(seq(1, n^2), n, n)
# i
ir <- cbind(i[ , -1], i[ , 1])
# ir
il <- cbind(i[ , n], i[ , -n])
# il
id <- rbind(i[-1, ], i[1, ])
# id
iu <- rbind(i[n, ], i[-n, ])
# iu
x <- matrix(1, n, n)
@
<<coding>>=
co <- (i + (1 - n %% 2) * col(i)) %% 2
# co
ico0 <- i[co == 0]
ico1 <- i[co == 1]
@

Now we want to do a parallel tempering setup.
For a start we need the whole unnormalized log density.
<<logh>>=
logh <- function(x, tau) sum((x * x[ir] * br + x * x[id] * bd) / tau)
@
Then we need to decide how how many ``helper'' chains we want.
We'll start with one.  And we need to decide what temperature
to use for the helper.
<<help-setup>>=
x.help <- x
tau.help <- 0.3
@

Define update.
<<block-gibbs>>=
block.gibbs <- function(x, tau) {
foo <- x[ir] * br + x[id] * bd + x[il] * br[il] + x[iu] * br[iu]
foo <- foo / tau
p <- 1 / (1 + exp(- 2 * foo))
x[ico0] <- as.numeric(runif(n^2 / 2) < p[ico0]) * 2 - 1
foo <- x[ir] * br + x[id] * bd + x[il] * br[il] + x[iu] * br[iu]
foo <- foo / tau
p <- 1 / (1 + exp(- 2 * foo))
x[ico1] <- as.numeric(runif(n^2 / 2) < p[ico1]) * 2 - 1
return(x)
}
@

Set random seed (if we always want the same results, otherwise omit).
<<seed>>=
set.seed(42)
@

\section{Run One}

<<do-one-setup>>=
nbatch <- 100
blen <- 1e3
accept <- 0
@
<<do-one>>=
batch <- matrix(NA, nbatch, 3 * n^2)
for (ibatch in 1:nbatch) {
    xbatch <- rep(0, 3 * n^2)
    for (iiter in 1:blen) {
        x <- block.gibbs(x, tau)
        x.help <- block.gibbs(x.help, tau.help)
        r <- logh(x.help, tau) + logh(x, tau.help) -
            logh(x, tau) - logh(x.help, tau.help)
        if (r > 0 || runif(1) < exp(r)) {
            foo <- x.help
            x.help <- x
            x <- foo
            accept <- accept + 1
        }
        xbatch <- xbatch + as.numeric(c(x, x * x[ir], x * x[id]))
    }
    batch[ibatch, ] <- xbatch / blen
}
accept <- accept / (nbatch * blen)
mu <- apply(batch, 2, mean)
mcse <- apply(batch, 2, sd) / sqrt(nbatch)
xmu <- matrix(mu[1:n^2], n, n)
xmcse <- matrix(mcse[1:n^2], n, n)
xxrmu <- matrix(mu[n^2 + 1:n^2], n, n)
xxrmcse <- matrix(mcse[n^2 + 1:n^2], n, n)
xxdmu <- matrix(mu[2 * n^2 + 1:n^2], n, n)
xxdmcse <- matrix(mcse[2 * n^2 + 1:n^2], n, n)
@
<<do-one-print>>=
accept
round(xmu, 3)
round(xmcse, 3)
round(xxrmu, 3)
round(xxrmcse, 3)
round(xxdmu, 3)
round(xxdmcse, 3)
@

We happened to luck out with the helper temperature, our first guess
giving an acceptance rate in the right range (at least the range specified
by the assignment, no theory is available that says what acceptance rate
is optimal).

The maximum
of all the MCSE is \Sexpr{round(max(xmcse, xxrmcse, xxdmcse), 4)},
so (assuming the MCSE were correct, which they are not)
we only need to run about
$\Sexpr{round(max(xmcse, xxrmcse, xxdmcse) / 0.001, 1)}^2$
times longer, that is
\verb@nbatch * blen@ equal to
<<bfoo>>=
ntodo <- nbatch * blen * (max(xmcse, xxrmcse, xxdmcse) / 0.001)^2
ntodo
ntodo.exp <- floor(log10(ntodo))
ntodo.frac <- round(ntodo / 10^ntodo.exp, 1)
@
\verb@nbatch * blen@ equal to
$\Sexpr{ntodo.frac} \times 10^{\Sexpr{ntodo.exp}}$,
divided in any way so the batches are long enough.

Examination of the autocorrelation plots for these shows that
\verb@blen@ should be at least 20 times longer.

\section{Run Two}

<<do-two>>=
blen <- 20 * blen
accept <- 0
<<do-one>>
<<do-one-print>>
@

The maximum
of all the MCSE is \Sexpr{round(max(xmcse, xxrmcse, xxdmcse), 4)},
so we only need to run about
$\Sexpr{round(max(xmcse, xxrmcse, xxdmcse) / 0.001, 1)}^2$
times longer, that is
\verb@nbatch * blen@ equal to
<<bfoo-too>>=
<<bfoo>>
@
\verb@nbatch * blen@ equal to
$\Sexpr{ntodo.frac} \times 10^{\Sexpr{ntodo.exp}}$,
divided in any way so the batches are long enough.

\section{Run Three}

<<do-three>>=
blen <- 2 * blen
accept <- 0
<<do-one>>
<<do-one-print>>
@

The maximum
of all the MCSE is \Sexpr{round(max(xmcse, xxrmcse, xxdmcse), 4)},
so we only need to run about
$\Sexpr{round(max(xmcse, xxrmcse, xxdmcse) / 0.001, 1)}^2$
times longer, that is
\verb@nbatch * blen@ equal to
<<bfoo-three>>=
<<bfoo>>
@
\verb@nbatch * blen@ equal to
$\Sexpr{ntodo.frac} \times 10^{\Sexpr{ntodo.exp}}$,
divided in any way so the batches are long enough.

\end{document}

\end{document}