%% ================================================================
\section{Introduction}

Coin flips are used as a motivating example to describe why one would want to
use the Metropolis-Hastings algorithm.  The algorithm is presented, illustrated
by example, and then proved correct.  See \cite{bib:kerl} for probability
terminology and notation used in this paper.

%% ================================================================
\section{Probability distributions}
\label{sec:four_models}

Consider flipping $n$ coins.  Encode a heads-up flip as $+1$ and a tails-up
flip as $-1$.  Then the result of the $i$th flip is $\omega_i \in \{+1,-1\}$.
I'll discuss four scenarios, with increasing levels of complexity.  Write
ordered $n$-tuples
$$
	\vecomega = (\omega_1, \ldots, \omega_n).
$$
Notice that there are $N=2^n$ outcomes in the probability space
$$
	\Omega = \{(\omega_1, \ldots, \omega_n) : \omega_i = \pm 1\}.
$$
For example, with $n=3$ coins, there are $N=8$ outcomes: HHH, HHT, HTH, HTT,
THH, THT, TTH, and TTT.

%% ----------------------------------------------------------------
\subsection{IID fair flips}
\label{sec:model_1}

In the first case I consider, the flips are independent and identically
distributed.  Furthermore, each coin is fair, which is to say
$$
	P_i(+1)=1/2 \qqand P_i(-1)=1/2.
$$
Due to independence, the probabilities factor:
\begin{align*}
	P(\vecomega) &=
		P(\omega_1, \ldots, \omega_n) \\
	&= P_1(\omega_1) \cdots P_n(\omega_n) \\
	&= \prod_{i=1}^n \frac{1}{2}
		= \left(\frac{1}{2}\right)^n.
\end{align*}
For example, with $n=3$, the probability of heads, tails, heads is
$$
	P(+1, -1, +1) = 1/8.
$$
The other 7 possible outcomes also occur with probability $1/8$.
Here's the PMF:

$$
\begin{tabular}{|rrrrr|c|c|}
\hline \multicolumn{5}{|c|}{($\omega_1,\omega_2,\omega_3$)}
	& $P(\omega_1,\omega_2,\omega_3)$ \\
\hline
\hline ( &  1 &  1 &  1 & ) & 0.125 \\
\hline ( &  1 &  1 & -1 & ) & 0.125 \\
\hline ( &  1 & -1 &  1 & ) & 0.125 \\
\hline ( &  1 & -1 & -1 & ) & 0.125 \\
\hline ( & -1 &  1 &  1 & ) & 0.125 \\
\hline ( & -1 &  1 & -1 & ) & 0.125 \\
\hline ( & -1 & -1 &  1 & ) & 0.125 \\
\hline ( & -1 & -1 & -1 & ) & 0.125 \\
\hline
\end{tabular}
$$

%% ----------------------------------------------------------------
\subsection{IID flips}
\label{sec:model_2}

In the second case I consider, the $n$ flips are still IID but the coin isn't
necessarily fair.  Each lands heads-up with probability $p$, where $0 \le p \le
1$.  Due to independence, the probabilities still factor.  Let
\begin{equation}
	\eta_i = \begin{cases}
		0, & \omega_i = +1 \\
		1, & \omega_i = -1.
	\end{cases}
\end{equation}
That is,
\begin{equation}
\label{eqn:eta}
	\omega_i = 1 - 2 \eta_i
	\qqand
	\eta_i   = \frac{1 - \omega_i}{2}.
\end{equation}
Then
\begin{align*}
	P(\vec\omega) &=
		P(\omega_1, \ldots, \omega_n) \\
	&= P_1(\omega_1) \cdots P_n(\omega_n) \\
	&= \prod_{i=1}^n p^{1-\eta_i} (1-p)^{\eta_i}.
\end{align*}
With $n=3$, the probability of heads, tails, heads is
\begin{align*}
	P(1, -1, 1) &=
	\Big(p^1(1-p)^0\Big) \;
	\Big(p^0(1-p)^1\Big) \;
	\Big(p^1(1-p)^0\Big) \\
	&= p^2(1-p).
\end{align*}

The PMF is as follows, with $p=0.9$:

$$
\begin{tabular}{|rrrrr|c|c|}
\hline \multicolumn{5}{|c|}{($\omega_1,\omega_2,\omega_3$)}
	& $P(\omega_1,\omega_2,\omega_3)$ \\
\hline
\hline ( &  1 &  1 &  1 & ) & 0.729 \\
\hline ( &  1 &  1 & -1 & ) & 0.081 \\
\hline ( &  1 & -1 &  1 & ) & 0.081 \\
\hline ( &  1 & -1 & -1 & ) & 0.009 \\
\hline ( & -1 &  1 &  1 & ) & 0.081 \\
\hline ( & -1 &  1 & -1 & ) & 0.009 \\
\hline ( & -1 & -1 &  1 & ) & 0.009 \\
\hline ( & -1 & -1 & -1 & ) & 0.001 \\
\hline
\end{tabular}
$$

%% ----------------------------------------------------------------
\subsection{Independent flips}
\label{sec:model_3}

In the third case I consider, the coin flips are still independent but they are
not identically distributed: each has its own probability $p_i$ of heads.  That
is,
$$
	P_i(\omega_i = +1) = p_i.
$$
Then
\begin{equation}
\label{eqn:P_eta}
\begin{aligned}
	P(\vecomega) &=
		P(\omega_1, \ldots, \omega_n) \\
	&= P_1(\omega_1) \cdots P_n(\omega_n) \\
	&= \prod_{i=1}^n p_i^{1-\eta_i} (1-p_i)^{\eta_i}.
\end{aligned}
\end{equation}
With $n=3$, the probability of heads, tails, heads is
\begin{align*}
	P(1, -1, 1) &= p_1 (1-p_2) p_3.
\end{align*}

With $p_1=0.6$, $p_2=0.7$, and $p_3=0.8$, the PMF looks like this:

$$
\begin{tabular}{|rrrrr|c|c|}
\hline \multicolumn{5}{|c|}{($\omega_1,\omega_2,\omega_3$)}
	& $P(\omega_1,\omega_2,\omega_3)$ \\
\hline
\hline ( &  1 &  1 &  1 & ) & 0.336 \\
\hline ( &  1 &  1 & -1 & ) & 0.084 \\
\hline ( &  1 & -1 &  1 & ) & 0.144 \\
\hline ( &  1 & -1 & -1 & ) & 0.036 \\
\hline ( & -1 &  1 &  1 & ) & 0.224 \\
\hline ( & -1 &  1 & -1 & ) & 0.056 \\
\hline ( & -1 & -1 &  1 & ) & 0.096 \\
\hline ( & -1 & -1 & -1 & ) & 0.024 \\
\hline
\end{tabular}
$$

%% ----------------------------------------------------------------
\subsection{Dependent flips}
\label{sec:model_4}

In the fourth and final case I consider, the coin flips are no longer
independent --- perhaps we are flipping magnetized coins, where the result of a
flip is influenced by one or more of its neighbors.  The PMF is a function with
potentially $2^n$ different input/output values.  There are all sorts of PMFs
one could write down to define a probability model for dependent coin flips.
For example:
$$
\begin{tabular}{|rrrrr|c|c|}
\hline \multicolumn{5}{|c|}{($\omega_1,\omega_2,\omega_3$)}
	& $P(\omega_1,\omega_2,\omega_3)$ \\
\hline
\hline ( &  1 &  1 &  1 & ) & 0.4 \\
\hline ( &  1 &  1 & -1 & ) & 0.0 \\
\hline ( &  1 & -1 &  1 & ) & 0.1 \\
\hline ( &  1 & -1 & -1 & ) & 0.1 \\
\hline ( & -1 &  1 &  1 & ) & 0.0 \\
\hline ( & -1 &  1 & -1 & ) & 0.1 \\
\hline ( & -1 & -1 &  1 & ) & 0.1 \\
\hline ( & -1 & -1 & -1 & ) & 0.2 \\
\hline
\end{tabular}
$$
However, it's not clear what kind of physical situation could lead to such an
ad-hoc PMF.

There is one particular family of models --- the Ising family --- which
contains enough flexibility to describe plenty of dependent-flip situations.
Furthermore, the Ising family has some physical plausibility.

%Here's one possible PMF.  (This is from a 1D Ising model as discussed in
%section \ref{sec:ising_1d}, with anti-aligning nearest-neighbor connections and
%a small downward magnetic field.)
%
%$$
%\begin{tabular}{|rrrrr|c|c|}
%\hline \multicolumn{5}{|c|}{($\omega_1,\omega_2,\omega_3$)}
%	& $P(\omega_1,\omega_2,\omega_3)$ \\
%\hline
%\hline ( &  1 &  1 &  1 & ) & 0.0038136 \\
%\hline ( &  1 &  1 & -1 & ) & 0.0420384 \\
%\hline ( &  1 & -1 &  1 & ) & 0.3106238 \\
%\hline ( &  1 & -1 & -1 & ) & 0.0627139 \\
%\hline ( & -1 &  1 &  1 & ) & 0.0420384 \\
%\hline ( & -1 &  1 & -1 & ) & 0.4633963 \\
%\hline ( & -1 & -1 &  1 & ) & 0.0627139 \\
%\hline ( & -1 & -1 & -1 & ) & 0.0126617 \\
%\hline
%\end{tabular}
%$$

%% ================================================================
\newpage
\section{The one-dimensional Ising model}
\label{sec:ising_1d}

I'll refer to \bfidx{spins} rather than \bfidx{flips}: heads and tails will be
replaced with spin-up and spin-down, respectively.  A sequence such as HHTH
will now be $++-+$, and I will call that a \bfidx{spin configuration}.

In the one-dimensional \plainidx{Ising model}, there are $n$ spins and $N=2^n$
possible states (spin configurations), and the individual spins in a
configuration are not necessarily independent.

Let $\vecomega=(\omega_1, \ldots, \omega_n)$ as before.  We define an
\bfidx{energy function} on the spins $\omega_i$ by
$$
	E(\vecomega) = \sum_{i=1}^n \sum_{j=1}^n S_{ij} \omega_i \omega_j
		+ \sum_{i=1}^n h_i \omega_i.
$$
The $S_{ij}$'s are \bfidx{coupling coefficients} which specify the
dependence of one spin on another.
The $h_i$'s are the \bfidx{magnetic field term}, which breaks the symmetry.

We get a probability model using the statistical-mechanical notation that
\emph{energy is the logarithm of probability}.  The probability of each
configuration $\vecomega$ is a \plainidx{Boltzmann distribution} with
\plainidx{inverse temperature} $\beta=1/kT$ (where $k$ is \plainidx{Boltzmann's
constant}):
\begin{equation}
\label{eqn:P_exp_beta_E}
	P(\vecomega) \propto e^{-\beta E(\vecomega)}.
\end{equation}
To normalize this, so that $\sum_{\vecomega} P(\vecomega)=1$,
we divide by
$$
	Z = {\sum_{k=1}^N e^{-\beta E(\vecomega^{(k)})}}
$$
where $\vecomega^{(k)}$ runs over all possible states.  This denominator has a
special name:  it is called a \bfidx{partition function}.  Thus,
$$
	P(\vecomega^{(j)}) = \frac{e^{-\beta E(\vecomega^{(j)})}}
		{\sum_{k=1}^N e^{-\beta E(\vecomega^{(k)})}}.
$$
Some handy notation: Write the \bfidx{probability numerator}
as
\begin{equation}
\label{eqn:P*}
	P^*(\vecomega^{(j)}) = e^{-\beta E(\vecomega^{(j)})}.
\end{equation}
Then
\begin{equation}
\label{eqn:PZ}
	P(\vecomega^{(j)}) = \frac{P^*(\vecomega^{(j)})}{Z}.
\end{equation}

We will see (if I get to it) that the diagonal terms $S_{ii}$ do not affect
spin probabilities.  Here are some examples of $S$ matrices, with
$n=3$.

$$
\begin{tabular}{ccc}
Non-interacting: &
Mean field: &
Aligning nearest-neighbor: \\
$
S =
\left(\begin{array}{rrr}
	0 & 0 & 0 \\
	0 & 0 & 0 \\
	0 & 0 & 0 \\
\end{array}\right)
$
&
$
S = 
\left(\begin{array}{rrr}
	0 & 1 & 1 \\
	1 & 0 & 1 \\
	1 & 1 & 0 \\
\end{array}\right)
$
&
$
S = 
\left(\begin{array}{rrr}
	0 & 1 & 0 \\
	0 & 0 & 1 \\
	0 & 0 & 0 \\
\end{array}\right)
$
\\
\\
\\
Anti-aligning nearest-neighbor: &
Aligning nearest-neighbor &
Anti-aligning nearest-neighbor \\
&
with periodic boundary conditions: &
with periodic boundary conditions: \\
$
S = 
\left(\begin{array}{rrr}
	0 & -1 &  0 \\
	0 &  0 & -1 \\
	0 &  0 &  0 \\
\end{array}\right)
$
&
$
S = 
\left(\begin{array}{rrr}
	0 & 1 & 0 \\
	0 & 0 & 1 \\
	1 & 0 & 0 \\
\end{array}\right)
$
&
$
S = 
\left(\begin{array}{rrr}
	0 & -1 &  0 \\
	0 &  0 & -1 \\
	-1&  0 &  0 \\
\end{array}\right)
$
\\
\end{tabular}
$$

Here are the PMFs for the first three of those, with $h=0$:
$$
\begin{tabular}{|rrrrr|c|c|c|}
\hline \multicolumn{5}{|c|}{($\omega_1,\omega_2,\omega_3$)}
	& $P(\omega_1,\omega_2,\omega_3)$
	& $P(\omega_1,\omega_2,\omega_3)$
	& $P(\omega_1,\omega_2,\omega_3)$  \\
\multicolumn{5}{|c|}{}
	& Non-interacting
	& Mean field
	& Aligning nearest neighbor \\
\hline
\hline ( &  1 &  1 &  1 & ) & 0.125 &  0.4994973 &  0.3879017 \\
\hline ( &  1 &  1 & -1 & ) & 0.125 &  0.0001676 &  0.0524968 \\
\hline ( &  1 & -1 &  1 & ) & 0.125 &  0.0001676 &  0.0071047 \\
\hline ( &  1 & -1 & -1 & ) & 0.125 &  0.0001676 &  0.0524968 \\
\hline ( & -1 &  1 &  1 & ) & 0.125 &  0.0001676 &  0.0524968 \\
\hline ( & -1 &  1 & -1 & ) & 0.125 &  0.0001676 &  0.0071047 \\
\hline ( & -1 & -1 &  1 & ) & 0.125 &  0.0001676 &  0.0524968 \\
\hline ( & -1 & -1 & -1 & ) & 0.125 &  0.4994973 &  0.3879017 \\
\hline
\end{tabular}
$$

%% ================================================================
\newpage
\section{Independence and CDF sampling}

In this section we discuss two sampling methods which are simpler than
Metropolis-Hastings:  independence sampling, which works for the independent
cases, and CDF sampling which works for all four examples of section
\ref{sec:four_models} but doesn't scale well as $n$ increases.

%% ----------------------------------------------------------------
\subsection{Sample proportions}
\label{sec:Q}

Suppose we flip four fair coins.  Getting all four heads, or all four tails, is
no huge surprise --- but getting two heads and two tails is more likely.  On
the other hand, if we flip a thousand fair coins, all heads or all tails are
very unlikely and we expect about 500 heads in general.  It is quite possible
to quantify these expectations precisely, in terms of sample mean and variance,
using the binomial theorem as described in \cite{bib:kerl}:  IID fair flips
lend themselves to neat theoretical analysis.  For purposes of this paper,
though, the goal is to describe experimental methods which are useful for
probability models which are \emph{not} theoretically neat.  It is more
important to look at sample proportions, which allow us to tabulate our
experimental results.

Suppose we conduct (in a manner which subsequent sections will make precise)
$M$ experiments of an $n$-coin flip, using any of the four models of section
\ref{sec:four_models}, with outcomes
$$
	\vecomega^{(1)}, \ldots, \vecomega^{(M)}.
$$
Then, for $j=1, \ldots, N=2^n$, define $Q_i$ to be the fraction of
$\vecomega$'s which landed in state $i$.  I also sometimes write, for brevity,
$$
	\bfP = (P_1, \ldots, P_N)
	\qand
	\bfQ = (Q_1, \ldots, Q_N).
$$

For example, if I flip $N=3$ fair coins $M=10$ times, I might get the
following:
$$
\begin{tabular}{|l|rrrrr|}
\hline $i$ & \multicolumn{5}{|c|}{$\vecomega^{(i)}$} \\
\hline
\hline  1 & ( &  1 &  1 &  1 & ) \\
\hline  2 & ( & -1 &  1 & -1 & ) \\
\hline  3 & ( & -1 &  1 &  1 & ) \\
\hline  4 & ( & -1 &  1 & -1 & ) \\
\hline  5 & ( & -1 & -1 &  1 & ) \\
\hline  6 & ( &  1 & -1 &  1 & ) \\
\hline  7 & ( &  1 & -1 &  1 & ) \\
\hline  8 & ( & -1 & -1 &  1 & ) \\
\hline  9 & ( &  1 & -1 & -1 & ) \\
\hline 10 & ( &  1 &  1 &  1 & ) \\
\hline
\end{tabular}
$$
Then I can tabulate the sample proportions, which I can compare against the
PMF:
$$
\begin{tabular}{|rrrrr|c|c|c|}
\hline \multicolumn{5}{|c|}{$\vecomega$}
	& $P(\vecomega)$
	& $Q(\vecomega)$ \\
\hline
\hline ( &  1 &  1 &  1 & ) & 0.125 & 0.2 \\
\hline ( &  1 &  1 & -1 & ) & 0.125 & 0.0 \\
\hline ( &  1 & -1 &  1 & ) & 0.125 & 0.2 \\
\hline ( &  1 & -1 & -1 & ) & 0.125 & 0.1 \\
\hline ( & -1 &  1 &  1 & ) & 0.125 & 0.1 \\
\hline ( & -1 &  1 & -1 & ) & 0.125 & 0.2 \\
\hline ( & -1 & -1 &  1 & ) & 0.125 & 0.2 \\
\hline ( & -1 & -1 & -1 & ) & 0.125 & 0.0 \\
\hline
\end{tabular}
$$
If I re-run this experiment with $M=10000$ trials, I might get the following
sample proportions:

$$
\begin{tabular}{|rrrrr|c|c|c|}
\hline \multicolumn{5}{|c|}{$\vecomega$}
	& $P(\vecomega)$
	& $Q(\vecomega)$ \\
\hline
\hline ( &  1 &  1 &  1 & ) & 0.125 & 0.1247 \\
\hline ( &  1 &  1 & -1 & ) & 0.125 & 0.1191 \\
\hline ( &  1 & -1 &  1 & ) & 0.125 & 0.1250 \\
\hline ( &  1 & -1 & -1 & ) & 0.125 & 0.1281 \\
\hline ( & -1 &  1 &  1 & ) & 0.125 & 0.1237 \\
\hline ( & -1 &  1 & -1 & ) & 0.125 & 0.1307 \\
\hline ( & -1 & -1 &  1 & ) & 0.125 & 0.1231 \\
\hline ( & -1 & -1 & -1 & ) & 0.125 & 0.1256 \\
\hline
\end{tabular}
$$

The \bfidx{law of large numbers} tells us that each of the eight $Q_i$'s
approach the respective $P_i$'s as the number $M$ of trials increases; the
\bfidx{central limit theorem} tells us about the distribution of the $Q_i$'s.

The next question is how to actually conduct such an experiment in a computer
program.

%% ----------------------------------------------------------------
\subsection{Sampling with the independence method}

For independent flips, e.g. any of the first three models in section
\ref{sec:four_models}, it's easy to produce samples.  Recall that $p_j$ is the
probability that the $j$th flip is heads, which we encode as $+1$.
Thus, the algorithm is straightforward.

\begin{tabbing}
\hspace{0.25in} \= \hspace{0.25in} \= \hspace{0.25in} \= \hspace{0.25in} \\
		\emph{Initialize the histogram bins:} \\
\+		For $i=1, \ldots, N$ (number of possible outcomes): \\
			counts${}_i = 0$. \\
\- \\
		\emph{Generate states and collect data:} \\
\+		For $i=1, \ldots, M$ (number of trials): \\
\+			For $j=1, \ldots, n$ (number of coins): \\
				$U=$ uniform random on $[0,1)$.  \\
\-				If $U < p_j$ then $\omega_j=+1$ else $\omega_j=-1$. \\
		Find the index $k$ of the state $\vecomega=(\omega_1, \ldots,
			\omega_n)$ which was just generated. \\
		Increment counts${}_k$ by 1. \\
\- \\
		\emph{Scale the histogram bins down to sample proportions:} \\
\+		For $i=1, \ldots, N$ (number of possible outcomes): \\
\-			Divide counts${}_i$ by $M$.
\end{tabbing}

The output of such an experiment was shown in section \ref{sec:Q}.

The main limitation of this algorithm is obvious:  it only works for
independent-flip models.

%% ----------------------------------------------------------------
\subsection{Sampling with the CDF method}

For all four models of section \ref{sec:four_models} --- independent or not ---
it is possible to sample using a technique which is almost as simple as the
independence sampling of the previous section.

The key concept is the \bfidx{cumulative distribution function}.  That is, we
write down the $N$ possible outcomes in some order $\vecomega^{(1)}$ through
$\vecomega^{(N)}$, then tabulate at the $j$th slot the probability of obtaining
states $\vecomega^{(1)}$ through $\vecomega^{(j)}$.  That is, the CDF, as a
function of $j$ from $1$ to $N$, is
$$
	C(j) = P(\{\vecomega^{(1)}, \ldots, \vecomega^{(j)}\}).
$$
E.g.
\begin{align*}
	C(1) &= P(\vecomega^{(1)} \\
	C(2) &= P(\vecomega^{(1)} +
		P(\vecomega^{(2)} \\
	C(3) &= P(\vecomega^{(1)} +
		P(\vecomega^{(2)} +
		P(\vecomega^{(3)} \\
	\vdots & \\
	C(N) &= P(\vecomega^{(1)} +
		\ldots +
		P(\vecomega^{(N)}.
\end{align*}

For example, using the IID-fair model, we have
$$
\begin{tabular}{|rrrrr|c|c|}
\hline \multicolumn{5}{|c|}{$\vecomega$}
	& $P(\vecomega)$
	& $CDF$ \\
\hline
\hline ( &  1 &  1 &  1 & ) & 0.125 & 0.125 \\
\hline ( &  1 &  1 & -1 & ) & 0.125 & 0.250 \\
\hline ( &  1 & -1 &  1 & ) & 0.125 & 0.375 \\
\hline ( &  1 & -1 & -1 & ) & 0.125 & 0.500 \\
\hline ( & -1 &  1 &  1 & ) & 0.125 & 0.625 \\
\hline ( & -1 &  1 & -1 & ) & 0.125 & 0.750 \\
\hline ( & -1 & -1 &  1 & ) & 0.125 & 0.875 \\
\hline ( & -1 & -1 & -1 & ) & 0.125 & 1.000 \\
\hline
\end{tabular}
$$
Likewise, using the aligning nearest-neighbor Ising model from section
\ref{sec:ising_1d}, with $h=0.1$, we have
$$
\begin{tabular}{|rrrrr|c|c|c|}
\hline \multicolumn{5}{|c|}{$\vecomega$}
	& $P(\vecomega)$
	& $CDF$ \\
\hline
\hline ( &  1 &  1 &  1 & ) &  0.3879017 & 0.3879017 \\
\hline ( &  1 &  1 & -1 & ) &  0.0524968 & 0.4403985 \\
\hline ( &  1 & -1 &  1 & ) &  0.0071047 & 0.4475032 \\
\hline ( &  1 & -1 & -1 & ) &  0.0524968 & 0.5000000 \\
\hline ( & -1 &  1 &  1 & ) &  0.0524968 & 0.5524968 \\
\hline ( & -1 &  1 & -1 & ) &  0.0071047 & 0.5596015 \\
\hline ( & -1 & -1 &  1 & ) &  0.0524968 & 0.6120983 \\
\hline ( & -1 & -1 & -1 & ) &  0.3879017 & 1.0000000 \\
\hline
\end{tabular}
$$

The algorithm is as follows.

\begin{tabbing}
\hspace{0.25in} \= \hspace{0.25in} \= \hspace{0.25in}
\= \hspace{0.25in} \= \hspace{0.25in} \\
%
	\emph{Turn the PMF into a CDF:} \\
	Sum = 0. \\
\+	For $i=1, \ldots, N$ (number of possible states): \\
		Sum := sum + $p_i$. \\
		CDF${}_i$ = sum. \\
%	num_states = len(Ps)
%	CDF = [0.0] * num_states
%	sum = 0.0
%	for k in range(0, num_states):
%		sum += Ps[k]
%		CDF[k] = sum
%	return CDF
\- \\
	\emph{Initialize the histogram bins:} \\
\+		For $i=1, \ldots, N$ (number of possible outcomes): \\
			counts${}_i = 0$. \\
\- \\
		\emph{Generate states and collect data:} \\
\+		For $i=1, \ldots, M$ (number of trials): \\
			$U=$ uniform random on $[0,1)$.  \\
\+			For $j=1, \ldots, N$: \\
\+				If $U < \textrm{CDF}_j$: \\
					Increment counts${}_j$ by 1; go the next iteration of the
					$i$ loop. \\
\-\-\- \\
	\emph{Scale the histogram bins down to sample proportions:} \\
\+		For $i=1, \ldots, N$ (number of possible outcomes): \\
\-			Divide counts${}_i$ by $M$.
\end{tabbing}

Example results are as follows, with $M=10000$ trials:
$$
\begin{tabular}{|rrrrr|r|r|r|}
\hline \multicolumn{5}{|c|}{$\vecomega$}
	& $P(\vecomega)$
	& $Q(\vecomega)$
	& Difference \\
\hline
\hline ( &  1 &  1 &  1 & ) & $0.3879017$ & $0.3832000$ & $ 0.0047017$ \\
\hline ( &  1 &  1 & -1 & ) & $0.0524968$ & $0.0480000$ & $ 0.0044968$ \\
\hline ( &  1 & -1 &  1 & ) & $0.0071047$ & $0.0070000$ & $ 0.0001047$ \\
\hline ( &  1 & -1 & -1 & ) & $0.0524968$ & $0.0506000$ & $ 0.0018968$ \\
\hline ( & -1 &  1 &  1 & ) & $0.0524968$ & $0.0557000$ & $-0.0032032$ \\
\hline ( & -1 &  1 & -1 & ) & $0.0071047$ & $0.0063000$ & $ 0.0008047$ \\
\hline ( & -1 & -1 &  1 & ) & $0.0524968$ & $0.0522000$ & $ 0.0002968$ \\
\hline ( & -1 & -1 & -1 & ) & $0.3879017$ & $0.3970000$ & $-0.0090983$ \\
\hline
\end{tabular}
$$

The limitation of this CDF-sampling algorithm is a bit subtle: it requires
computing probabilities for all $N$ states.  In the coin-flip example, $N=2^n$.
For, say, $n=100$, $2^n \approx 10^{69}$ which is (as I recall) on the order of
the number of protons in the universe.  One cannot even write down all possible
states.  Yet 100 coins isn't a lot so it seems like we should be able to do
\emph{something}.

There is even worse news.  Computing the \plainidx{partition function} $Z$
requires a sum over all $N$ states, which is also infeasible.  Not only can we
not list the probability of all states, we can't even compute the (normalized)
probability of even \emph{one} state, since (from equations \ref{eqn:P*} and
\ref{eqn:PZ}) $P(\vecomega)=P^*(\vecomega)/Z$.  It would be nice if we had
something which only depended on the probability numerator $P^*(\vecomega) =
e^{-\beta E(\vecomega)}$, which is easy to compute.

We will see in the following sections that the Metropolis-Hastings algorithm
overcomes both these problems.  First we will see what that algorithm is and
what it does; then, we'll see why it works.

%% ================================================================
\newpage
\section{Sampling with Metropolis-Hastings}
\label{sec:mh_sampling}

Here we present the essential components of the Metropolis-Hastings algorithm,
in pseudocode and a worked example.  Important, non-negligible practical
considerations are deferred until section \ref{sec:mh_choices}; correctness is
proved in section \ref{sec:MCMC}.

%% ----------------------------------------------------------------
\subsection{The algorithm}
\label{sec:mhalgo}

The goal is to select $M$ samples $\vecomega^{(1)}, \ldots, \vecomega^{(M)}$
from the PMF of the 1D Ising model with specified $S$, $h$, and $\beta$. The
states have probability numerator $P^*(\vecomega) = \exp(-\beta E(\vecomega))$
as given by equation \ref{eqn:P*}.  Here, I display spins as up (filled) or
down (hollow):

\begin{center}
\psfragscanon
\includegraphics[scale=1.0]{figures/ising1.eps}
\end{center}

\begin{itemize}

\item There is an \emphidx{inverse temperature} $\beta$:  For physical systems,
this has some meaning and, presumably, would be varied under some controlled
pattern.  For the independent models, $\beta$ cancels and has no effect (as
we'll see in section \ref{sec:DE_indep}) so you can give it any non-zero
value, e.g.  $\beta=1$.

\item Pick an \emph{initial configuration}.  Typically, there are three
choices: (1) Start with all spins down, i.e. $\vecomega = (-1, \ldots, -1)$.
(2) Start with all spins up, i.e. $\vecomega = (+1, \ldots, +1)$.  (3) Start
with $\vecomega$ selected from a uniform probability distribution on $\Omega$.
(See section \ref{sec:pup} for more information on choice of initial state.)

\item Select a \emph{site} $j$ and decide whether to flip $\omega_j$ to
$-\omega_j$.  (See section \ref{sec:rmk} for more information on site
selection.)

\begin{center}
\psfragscanon
\includegraphics[scale=1.0]{figures/ising2.eps}
\end{center}

\item This decision is made using the \emph{Metropolis prescription},
namely:

	\begin{itemize}
	\item Compute the change in energy $\Delta E = E(\vecomega')
		- E(\vecomega)$ which would be obtained if $\vecomega$ were
		sent to $\vecomega'$ by flipping $\omega_j$.
	\item One may compute $\Delta E$ by separately computing
		$E(\vecomega')$ and $E(\vecomega)$ and subtracting the two.  However,
		since the only change is at the site $j$, one may do some ad-hoc
		algebra (sections \ref{sec:DE_indep} and \ref{sec:DE_ising}) to derive
		an expression for $\Delta E$ which is less computationally expensive.
	\item Accept the change with probability
		$$
			\min\{1, e^{-\Delta E}\}.
		$$
		Otherwise, the change is said to be rejected.
	\end{itemize}

	This is called a \emph{Metropolis step}\index{step}.

\item Looping through all $n$ sites from $j=1$ to $j=n$, performing a
Metropolis step at each site $j$, is called a \emph{Metropolis
sweep}\index{sweep}.

\item Doing $M$ sweeps is a complete experiment.  (See section
\ref{sec:M} for more information on choice of $M$.)

\item If one realizes a random variable $X(\vecomega)$ at each of $M$ sweeps,
averaging $X$ over the $M$ sweeps, one obtains an approximation $\overline X$
for the expectation $E[X]$.  (See section \ref{sec:1D_Ising_RVs} for more
information on choice of random variable.)

\end{itemize}

The above may be rephrased in pseudocode as follows:

\begin{tabbing}
\hspace{0.25in} \= \hspace{0.25in} \= \hspace{0.25in} \= \hspace{0.25in} \\
		\emph{Pick the initial state $\vecomega$:} \\
\+		For $i=1, \ldots, n$ (number of spins): \\
\-			$\omega_i = 0$ (spins-down) \\
		OR\\
\+		For $i=1, \ldots, n$ (number of spins): \\
\-			$\omega_i = 1$ (spins-up) \\
 		OR\\
\+		For $i=1, \ldots, n$ (number of spins): \\
			$\omega_i = \pm 1$ with probabilities $1/2,1/2$
			(uniform random initial state). \\
\- \\
		\emph{Generate states and collect data:} \\
\+		For $i=1, \ldots, M$ (number of trials): \\
			\emph{One Metropolis sweep:} \\
\+			For $j=1, \ldots, n$ (number of spins): \\
				\emph{One Metropolis step:} \\
				$\Delta E = E(\vecomega')-E(\vecomega)$ 
					where $\vecomega'$ has opposite spin at site $j$. \\
				Let $V=\min\{1,e^{-\beta \Delta E}\}$. \\
				Generate $U$ uniformly distributed on $[0,1)$. \\
\+				If $U < V$: \\
\-					Accept the new state: $\vecomega := \vecomega'$. \\
\+				Else: \\
\-					Reject the new state:  keep $\vecomega$. \\
		Find the index $k$ of the state $\vecomega=(\omega_1, \ldots,
			\omega_n)$ which was just generated. \\
		Increment counts${}_k$ by 1.  (Do this whether the new state was
		accepted or not.) \\
\-\- \\
		\emph{Scale the histogram bins down to sample proportions:} \\
\+		For $k=1, \ldots, N$ (number of possible outcomes): \\
\-			Divide counts${}_k$ by $M$.
\end{tabbing}

Notice that a \emph{step} of this algorithm touches one of the $n$ sites; a
\emph{sweep} touches all $n$ sites.  We collect data only once per sweep, since
only then have all sites been given a chance to change.

This is the basic algorithm.  Important caveats, though, are presented in
section \ref{sec:mh_choices}.

%% ----------------------------------------------------------------
\subsection{Walk-through}

Here are two sweeps with $n=3$, $\beta=1$, aligning nearest-neighbor model,
initial state $(1,1,1)$, and $h=0.1$.  These are the energies and probability
numerators:

$$
\begin{tabular}{|c|rrrrr|c|c|c|}
\hline State number & \multicolumn{5}{|c|}{$\vecomega$}
	& $E(\vecomega)$
	& $P^*(\vecomega)$ \\
\hline
\hline 1 & ( &  1 &  1 &  1 & ) &   -1.7 & 5.4739474 \\
\hline 2 & ( &  1 &  1 & -1 & ) &    0.1 & 0.9048374 \\
\hline 3 & ( &  1 & -1 &  1 & ) &    2.1 & 0.1224564 \\
\hline 4 & ( &  1 & -1 & -1 & ) &   -0.1 & 1.1051709 \\
\hline 5 & ( & -1 &  1 &  1 & ) &    0.1 & 0.9048374 \\
\hline 6 & ( & -1 &  1 & -1 & ) &    1.9 & 0.1495686 \\
\hline 7 & ( & -1 & -1 &  1 & ) &   -0.1 & 1.1051709 \\
\hline 8 & ( & -1 & -1 & -1 & ) &   -2.3 & 9.9741825 \\
\hline
\end{tabular}
$$

\begin{verbatim}
Initial state:       ( 1  1  1 )

Sweep 1:
    omega       j  state#  Delta E  U         V         Accept/reject New E
  ------------ --- ------  -------  --------- --------- ------------- -----
  ( 1  1  1 )   1   1      1.8      0.9731406 0.1652989 Reject        -1.7
  ( 1  1  1 )   2   1      3.8      0.3204335 0.0223708 Reject        -1.7
  ( 1  1 -1 )   3   2      1.8      0.1590330 0.1652989 Accept         0.1
State after sweep 1: ( 1  1 -1 )

Sweep 2:
    omega       j  state#  Delta E  U         V         Accept/reject New E
  ------------ --- ------  -------  --------- --------- ------------- -----
  ( 1  1 -1 )   1   2      1.8      0.6781967 0.1652989 Reject         0.1
  ( 1 -1 -1 )   2   4     -0.2      0.9644570 1.0000000 Accept        -0.1
  ( 1 -1 -1 )   3   4      2.2      0.2357908 0.1108032 Reject        -0.1
State after sweep 2: ( 1 -1 -1 )
\end{verbatim}

%% ----------------------------------------------------------------
\subsection{Numerical results}

Example results are as follows, $n=3$, aligning nearest-neighbor model, $h=0.1$,
and $M=10000$ trials:
$$
\begin{tabular}{|rrrrr|r|r|r|}
\hline \multicolumn{5}{|c|}{$\vecomega$}
	& $P(\vecomega)$
	& $Q(\vecomega)$
	& Difference \\
\hline
\hline ( &  1 &  1 &  1 & ) &   $0.3879017$ & $0.3875000$ & $ 0.0004017$ \\
\hline ( &  1 &  1 & -1 & ) &   $0.0524968$ & $0.0514000$ & $ 0.0010968$ \\
\hline ( &  1 & -1 &  1 & ) &   $0.0071047$ & $0.0066000$ & $ 0.0005047$ \\
\hline ( &  1 & -1 & -1 & ) &   $0.0524968$ & $0.0523000$ & $ 0.0001968$ \\
\hline ( & -1 &  1 &  1 & ) &   $0.0524968$ & $0.0530000$ & $-0.0005032$ \\
\hline ( & -1 &  1 & -1 & ) &   $0.0071047$ & $0.0070000$ & $ 0.0001047$ \\
\hline ( & -1 & -1 &  1 & ) &   $0.0524968$ & $0.0523000$ & $ 0.0001968$ \\
\hline ( & -1 & -1 & -1 & ) &   $0.3879017$ & $0.3899000$ & $-0.0019983$ \\
\hline
\end{tabular}
$$
% The initial state was spins-up, sweeps were sequential, and thermalization was
% declared after 10000 pre-trials.

%section \ref{sec:mh_choices}

% xxx xref fwd ...

%% ----------------------------------------------------------------
\subsection{Heuristics}

Recall from section \ref{sec:mhalgo} that the Metropolis \plainidx{transition
probability} is
$$
	\min\{1.0, \exp(-\beta \Delta_E)\}.
$$
So, the Metropolis algorithm certainly transitions to states with non-negative
$\exp(-\beta \Delta E)$, and randomly transitions to states with negative
$\exp(-\beta \Delta E)$.  The certain transition occurs when
\begin{align*}
	\exp(-\beta \Delta E) &\ge 1 \\
	-\beta \Delta E &\ge \log(1) = 0 \\
	\beta \Delta E &\le 0 \\
	\Delta E &\le 0.
\end{align*}
Thus, we may also say that the Metropolis algorithm certainly transitions to
lower-energy or same-energy states, and randomly transitions to higher-energy
states.

$$
\begin{tabular}{|c|c|c|}
\hline && \\
Certain transition & $\exp(-\beta \Delta E) \ge 1$ & $\Delta E \le 0$ \\
&& \\
\hline && \\
Random transition & $\exp(-\beta \Delta E)  <  1$ & $\Delta E  >  0$ \\
&& \\
\hline
\end{tabular}
$$

%% ----------------------------------------------------------------
\subsection{Random variables for 1D Ising}
\label{sec:1D_Ising_RVs}

Here, the RVs are counts $Q_i$ as described in section \ref{sec:Q}.  This is
just for didactic purposes.  Often people look at \emph{sample means}
(averages) of some quantity, taken at each sample.  One can also define energy,
total magnetization, correlations between specified pairs of sites, etc.

%% ================================================================
\newpage
\section{Design choices and caveats for Metropolis-Hastings}
\label{sec:mh_choices}

For clarity, I presented the basic algorithm in section \ref{sec:mhalgo}.  Here
are additional practical considerations which must be addressed in any serious
use of Metropolis-Hastings.  Correctness of the algorithm is proved in section
\ref{sec:MCMC}.

%% ----------------------------------------------------------------
\subsection{Number of trials}
\label{sec:M}

Choice of the number of iterations $M$:  How to do this wisely is a story
in itself.  For this paper, I suggest trying a few runs with each $M$, for
varying $M$'s.  For example, do half a dozen runs with $M=10000$, another
half-dozen runs with $M=20000$, etc.  Then choose an $M$ which reduces the
standard deviation of the sample mean to your desired number of decimal places.

%% ----------------------------------------------------------------
\subsection{Initial state}
\label{sec:pup}

Choice of the initial state is as described in section \ref{sec:mhalgo}:
typically lowest-energy/highest-energy (e.g. here the all-down/all-up states),
or random configuration.  Results should be independent of initial state.
Failure to achieve this may be due to inadequate thermalization (section
\ref{sec:thermo}).

%% ----------------------------------------------------------------
\subsection{Sequential vs. random sweeps}
\label{sec:rmk}

The decision here is whether to sweep through all sites from (say) left to
right, or to move through them randomly.  There is some debate in the
literature on this topic (for which I lack citations here).

%% ----------------------------------------------------------------
\subsection{Thermalization}
\label{sec:thermo}

One should first run $L$ Metropolis sweeps of the system, discarding any data
collected, before running the $M$ sweeps in which data are accumulated.  The
$L$ sweeps are called the \emph{thermalization phase}; the $M$ sweeps are
called the \emph{accumulation phase}.

xxx clarify with numerical/graphical examples (describe in terms of $P$ plots
in section xxx):  The underlying concern is the convergence of the Metropolis
probability distribution to the stable distribution of its implicit Markov
chain.

Put up $E$ plots?

According to Tom Kennedy, there is no general method to determine whether the
system has thermalized.  xxx discuss $E$-plot idea (if it works out ...);
manual approaches.

%% ----------------------------------------------------------------
\subsection{Separation of states}

xxx highly annp --- the distribution is supported on all-spins-up and
all-spins-down.  The Metro one-at-a-time method will always give rejects.
Generalize, and cite some further reading.

%% ----------------------------------------------------------------
\subsection{Downsampling to avoid correlation}

Downsampling and batched means are discussed thoroughly in appendix B of
\cite{bib:kerldis}.  The discussion there also includes information about
Berg's autocorrelation detection.

%% ----------------------------------------------------------------
\subsection{$\Delta E$ computations for the independent model}
\label{sec:DE_indep}

As noted in section \ref{sec:mhalgo}, it is usually possible to compute $\Delta
E$ at less computational cost than $E(\vecomega')-E(\vecomega)$.  That is,
computing $E(\vecomega)$ requires visiting all $n$ sites of a spin
configuration, as does computing $E(\vecomega')$, but since the Metropolis
algorithm only modifies \emph{one} spin, we can skip most of the redundant
values which will only cancel out anyway.

When you use a Metropolis-Hastings on a new problem, you will likely be
doing some calculations of the type shown here.

xxx note $\beta$ shouldn't affect probabilities --- it shouldn't matter whether
we're flipping hot coins or cold ones.  then do the algebra and show why this
intuition is right.

xxx independence from temperature.

Using equations \ref{eqn:P_exp_beta_E} and \ref{eqn:P_eta},
\begin{align*}
	P(\vecomega) &\propto e^{-\beta E(\vecomega)} \\
	xxx & \\
	E(\vecomega) &= -\frac{1}{\beta} \log(P(\vecomega)) \\
	&= -\frac{1}{\beta} \log\left(
		\prod_{i=1}^n p_i^{1-\eta_i} (1-p_i)^{\eta_i}
	\right) \\
%
	&= -\frac{1}{\beta}
	\sum_{i=1}^n \left( (1-\eta_i)\log(p_i) + \eta_i\log(1-p_i) \right).
\end{align*}
xxx flip only at $j$th slot:  $\omega'_j=-\omega_j$ and so $\eta'_j=1-\eta_j$:
\begin{align*}
	\Delta E
	&= -\frac{1}{\beta} \left(
		\eta_j\log(p_j) + (1-\eta_j)\log(1-p_j)
		-(1-\eta_j)\log(p_j) - \eta_j\log(1-p_j)
	\right) \\
%
	&= -\frac{1}{\beta} \left(
		(1-2\eta_j)\log(1-p_j) -(1-2\eta_j)\log(p_j)
	\right) \\
%
	&= -\frac{1}{\beta}
		(1-2\eta_j)\log\left(\frac{1-p_j}{p_j}\right) \\
%
	&= \frac{\omega_j}{\beta} \log\left(\frac{p_j}{1-p_j}\right).
\end{align*}

xxx (or perhaps just state that the $\beta$'s cancel):
\begin{align*}
	\exp\left( -\beta \Delta E \right)
	&= \exp\left( - \omega_j \log\left(\frac{p_j}{1-p_j}\right) \right)
%
	= \left(\frac{1-p_j}{p_j}\right)^{\omega_j}.
\end{align*}

%xxx certain-transition case is
%\begin{align*}
%	\left(\frac{1-p_j}{p_j}\right)^{\omega_j} &\ge 1.
%\end{align*}
%
%xxx case $\omega_j=1$:  certain-transition case is
%\begin{align*}
%	\frac{1-p_j}{p_j} &\ge 1 \\
%	1-p_j &\ge p_j \\
%	1 &\ge 2p_j \\
%	1/2 &\ge p_j \\
%	p_j &\le 1/2.
%\end{align*}
%
%xxx case $\omega_j=-1$:  certain-transition case is
%\begin{align*}
%	\frac{p_j}{1-p_j} &\ge 1 \\
%	p_j &\ge 1-p_j \\
%	2p_j &\ge 1 \\
%	p_j &\ge 1/2.
%\end{align*}

xxx certain-transition case is
\begin{align*}
	\left(\frac{1-p_j}{p_j}\right)^{\omega_j} &\ge 1.
\end{align*}
Thus, when $\omega_j=1$, there is a certain transition when $p_j\le 1/2$; when
$\omega_j=1$, there is a certain transition when $p_j\ge 1/2$.  Since $p$ was
$P(\omega_j=1)$ --- the probability of heads at the $j$th slot --- the
Metropolis algorithm certainly sends individual spins to their more likely
states, and randomly sends them to their less likely states.

%% ----------------------------------------------------------------
\subsection{$\Delta E$ computations for the 1D Ising model}
\label{sec:DE_ising}

As discussed at the beginning of the previous section, and as mentioned in
section \ref{sec:mhalgo}, we can compute $\Delta E$ at less computational cost
than $E(\vecomega')-E(\vecomega)$.

The current state is $\vecomega$ and the energy is
$$
	E = -\sum_{ij} S_{ij} \omega_i \omega_j + \sum_i h_i \omega_i.
$$
The candidate new state $\vecomega'$ has opposite spin $\omega_k$.  For example,
with $n=4$,
$k=2$, and
$$
	\vecomega = (+1,+1,+1,+1),
$$
we have (numbering states from $0$ up to $n-1$):
$$
	\vecomega' = (+1,+1,-1,+1).
$$
Then the new energy is
$$
	E' = -\sum_{ij} S_{ij} \omega_i' \omega'_j + \sum_i h_i \omega_i'
$$
and the energy change is
\begin{align*}
	\Delta E &= E' - E \\
	&=	-\sum_{ij} S_{ij} \omega'_i \omega'_j
		+\sum_{ij} S_{ij} \omega_i \omega_j
		+\sum_i h_i (\omega'_i-\omega_i).
\end{align*}
Now, in my computer program I'm going to be computing this energy difference
very many times, so I want to omit as much redundant computation as possible.
For $i\ne k$ we have $\omega'_i = \omega_i$, and for $i=k$ we have
$\omega'_k=-\omega_k$.  So the last term is
$$
	h_k (\omega'_k-\omega_k) =
	-2 h_k \omega_k.
$$
The sums split into four pieces:
(1) $i\ne k$ and $j\ne k$,
(2) $i=k$ and $j\ne k$,
(3) $i\ne k$ and $j=k$, and
(4) $i=k$ and $j=k$.

The first piece is
$$
	-\sum_{i,j\ne k} S_{ij} \omega'_i \omega'_j
	+\sum_{i,j\ne k} S_{ij} \omega_i \omega_j
	=
	-\sum_{i,j\ne k} S_{ij} \omega_i \omega_j
	+\sum_{i,j\ne k} S_{ij} \omega_i \omega_j
	= 0.
$$

The second piece is
$$
	-\sum_{j\ne k} S_{kj} \omega'_k \omega'_j
	+\sum_{j\ne k} S_{kj} \omega_k \omega_j
	=
	 \sum_{j\ne k} S_{kj} \omega_k \omega_j
	+\sum_{j\ne k} S_{kj} \omega_k \omega_j
	=
	2\sum_{j\ne k} S_{kj} \omega_k \omega_j.
$$

The third piece is
$$
	-\sum_{i\ne k} S_{ik} \omega'_i \omega'_k
	+\sum_{i\ne k} S_{ik} \omega_i \omega_k
	=
	2\sum_{i\ne k} S_{ik} \omega_i \omega_k.
$$

The fourth piece is
$$
	-S_{kk} \omega'_k \omega'_k
	+S_{kk} \omega_k \omega_k
	=
	-S_{kk} \omega_k \omega_k
	+S_{kk} \omega_k \omega_k
	= 0.
$$
Thus the energy difference is
$$
	 2\left(\sum_{j\ne k} S_{kj} \omega_k \omega_j
	+\sum_{i\ne k} S_{ik} \omega_i \omega_k
	- h_k \omega_k\right).
$$

%% ================================================================
\newpage
\section{Markov chains}
\label{sec:markov_chains}

In order to prove the correctness of Metropolis-Hastings, we first need to
understand Markov chains.  We get \plainidx{Markov chains} by putting together
two concepts: \bfidx{stochastic processes} and \bfidx{conditional probability}.
I preface discussion of Markov chains with the \bfidx{die-tipping example}
which helps motivate the discussion.

Throughout, I assume elementary probability terminology as in
\cite{bib:kerl,bib:lawler,bib:gs}.  Namely: \plainidx{sample space};
\plainidx{probability measure}; \plainidx{probability mass function}
(\plainidx{PMF}); \plainidx{probability density function} (\plainidx{PDF});
\plainidx{random variable}.

%% ----------------------------------------------------------------
\subsection{The die-tipping experiments}
\label{sec:die_tip}

Most of this document uses coin-flipping and the 1D Ising model as a running
example.  Here, though, I describe two even simpler experiments (experiment A
and experiment B), which involve tipping a single die.  As simple as these are,
they encapsulate all the central points I want to make about \bfidx{Markov
chains}.  I'll return to coin flips and the Ising model in section xxx.

First recall that the pips on opposite faces of a die add to 7:
\begin{center}
\includegraphics[scale=0.3]{figures/single_die.eps}
\end{center}

\textbf{Experiment A}:  I set the die on the table with the one-face up.  Then
I close my eyes momentarily and tip it, so the one-face is now on the side.  I
then close my eyes and tip it again, and so on, until I've tipped it 20 times.

If I repeat this experiment --- placing the side and doing the 20 tips --- many
times, what's the probability distribution for the initial state, before
tipping?  I always start with 1 on top, so we have
$$
	\vecP^{(0)} = (1,0,0,0,0,0).
$$
(The six-tuple notation is the probability of 1-face up, 2-face up, and so on.
The superscript on the $\vecP$ indicates how many tips have occurred.)

What are the possibilities for the new top face after the first tip?  The
1-face was tipped over, and the 6-face was on the bottom so it can't be on top.
But 2, 3, 4, and 5 are equally likely.  So the probability distribution is
$$
	\vecP^{(1)} =
	\left(0,\frac{1}{4},\frac{1}{4},\frac{1}{4},\frac{1}{4},0\right).
$$
Note in particular that the up-face after each tip is not independent of the
up-face before that tip:  \bfidx{given} that the current up-face is, say, 3,
there are \bfidx{transition probabilities} for what the next up-face will be.

What about $\vecP^{(2)}$, $\vecP^{(3)}$, and so on?  It's an elementary, but
somewhat tedious, exercise in cases to show that
$$
	\vecP^{(2)} = \left(
	\frac{1}{4},\frac{1}{8},\frac{1}{8},\frac{1}{8},\frac{1}{8},\frac{1}{4}
	\right).
$$
One can continue this pattern, but one of my goals is to make an advertisement
for the computational advantages of Markov matrices.  We'll find out in section
\ref{sec:die_tip_2} how to easily compute $\vecP^{(k)}$.

Even before computing that, though, I certainly expect that after many tips,
the die won't remember that the starting position was with the 1-face up.  So,
$\vecP^{(k)}$ ought to begin to look uniform.

\textbf{Experiment B}:  The only difference here is that instead of starting
with the 1-face up, I instead first roll the die.  Then the initial
distribution is
$$
	\vecP^{(0)} = \left(
	\frac{1}{6},\frac{1}{6},\frac{1}{6},\frac{1}{6},\frac{1}{6},\frac{1}{6}
	\right).
$$
Again, you can go through the cases and find that
$$
	\vecP^{(1)} = \left(
	\frac{1}{6},\frac{1}{6},\frac{1}{6},\frac{1}{6},\frac{1}{6},\frac{1}{6}
	\right).
$$
and so on.  Again, Markov matrices in section \ref{sec:die_tip_2} will make this easier.

In experiment A, the probability distribution started out supported only on
1-up but (we expect) begins to approach uniform; in experiment B, the
probability distribution is already uniform and stays that way.

Key points:
\begin{itemize}
\item We have a choice of \bfidx{initial state} for the die-tipping
	experiment:  for example, one-up, or roll first, or others we could
	have used.
\item We have a \bfidx{probability distribution} $\vecP^{(0)}$ for the
	initial state.
\item We have probability distributions $\vecP^{(k)}$ for the
	\bfidx{subsequent states}.  These distributions evolve in $k$.
\item We have \bfidx{transition probabilities} constraining which face can
	appear on each step of the experiment.  These connect the probability
	distributions $\vecP^{(k)}$ and $\vecP^{(k+1)}$.  In section \ref{sec:die_tip_2} we will
	describe the evolution and the transition probabilities in terms of
	\bfidx{Markov chains}.
\item Some distributions evolve, or \plainidx{converge}, toward a
	\bfidx{limiting distribution}; this is called \bfidx{mixing} or
	\bfidx{thermalization}.  The limiting distribution (here, uniform) doesn't
	evolve.  This limiting distribution is called a \bfidx{stationary
	distribution} for the experiment.
\item For this experiment, it looks various distributions evolve toward the
	uniform distribution.  We'll see in section xxx how, and when, this happens.
\end{itemize}

%% ----------------------------------------------------------------
\subsection{Stochastic processes}

\begin{defn} Given a probability space $\Omega$ and a sequence $P^{(k)}$ of
probability measures for $k=0,1,2,\ldots$, a \bfidx{discrete-time stochastic
process}\index{stochastic process} is a sequence of random
variables\footnote{For the measure-theoretically inclined reader, $X^{(k)}$ is
a measurable function from $(\Omega,\mcF,P^{(k)})$ to $(S,\mcG)$ where $\mcF$
and $\mcG$ are sigma algebras on $\Omega$ and $S$, respectively.} $X^{(k)}$
from $\Omega$ to a \bfidx{state space} $S$.  \end{defn}

For concreteness, think of $\Omega$ as all $N=2^n$ heads/tails configuration of
$n$ coin flips, or all $N=6$ results of the landing of a single die, where we
allow the probability distributions $P^{(k)}$ to vary in time.

The state space might be another copy of $\Omega$, or $\R$, etc.  For example,
$X^{(k)}:\Omega \to \R$ might count the number of heads in a set of $n$ coin
flips.  Or, $W^{(k)}:\Omega \to \Omega$ might be the identity function.

(I should point out that many, if not most, authors define a stochastic process
with one ``long $\Omega$'', e.g. the probablity space is all \emph{sequences}
of $n$-coin-flip experiments.)

We categorize stochastic processes by their inputs and outputs:
\begin{itemize}

\item If $k=0,1,2,\ldots$, as above, then we have a \bfidx{discrete-time
stochastic process}, consisting of countably many random variables $X^{(k)}$.
If $k\in[0,+\infty)$, then we have a \bfidx{continuous-time stochastic
process}, consisting of uncountably many random variables $X^{(k)}$.

\item If each $X^{(k)}$ has only finitely or countably infinitely many possible
output values, then we have a \bfidx{discrete-valued stochastic process}, and
the distribution of $X^{(k)}$ is a PMF.  If each $X^{(k)}$ is a continuous
random variable, then we have a \bfidx{continuous-valued stochastic process},
and the distribution of $X^{(k)}$ is a PDF.

\end{itemize}

For this paper, I am considering discrete-time, discrete-valued stochastic
processes; Brownian motion is an example of a continuous-time,
continuous-valued stochastic process.

%% ----------------------------------------------------------------
\subsection{Conditional probability}

Forgetting for just a moment about stochastic process, consider two discrete
random variables $X$ and $Y$ with \plainidx{joint PMF}
\begin{align*}
	f_{X,Y}(x,y) &= P(X=x,Y=y).
\end{align*}
We can sum over $x$'s or $y$'s to obtain the \plainidx{marginal distributions}
of $X$ and $Y$:
\begin{align*}
	P(X=x_i) &= \sum_j P(Y=y_j,X=x_i)
	\qqand
	P(Y=y_j) = \sum_i P(Y=y_j,X=x_i).
\end{align*}
Recall the definition of \plainidx{conditional probability}:
\begin{align*}
	P(Y=y_j) &= \frac{P(Y=y_j,X=x_i)}{P(X=x_i)}
		\quad\textrm{(if $P(X=x_i)\ne 0$, otherwise 0)}.
\end{align*}
Multiply both sides by the denominator and sum over $i$:
\begin{align*}
	P(Y=y_j,X=x_i) &= P(Y=y_j\mid X=x_i) P(X=x_i).
\end{align*}
Then we have
\begin{align}
\label{eqn:cond_prob}
	P(Y=y_j) &= \sum_i P(Y=y_j\mid X=x_i) P(X=x_i).
\end{align}

For example, take $X$ to be the starting up-face of the die in experiment A or
B of section \ref{sec:die_tip}, and take $Y$ to be the up-face after the first
tip.  Then the intuition is that equation \ref{eqn:cond_prob} is the following.
How likely is it that $Y$ will be in state $j$?  Add up all the likelihoods of
all previous states $i$, weighted by the\index{transition probability}
probability that state $i$ transitions to state $j$.

In fact, if $X$ and $Y$ take finitely many values, say $n$ each, we can write
down an even nicer equation.  Take $n=2$ for example.  Then equation
\ref{eqn:cond_prob}, for $j=1,2$ simultaneously, is just a matrix product:
\begin{align*}
(P(Y=y_1), P(Y=y_2))
&=
(P(X=x_1), P(X=x_2))
\pcmatrix {
	P(Y=y_1 \mid X=x_1) & P(Y=y_1 \mid X=x_1) \\
	P(Y=y_1 \mid X=x_1) & P(Y=y_1 \mid X=x_1) \\
}
\end{align*}
For shorthand, write
$$
	P_i = P(X=x_i),
	\qquad
	P'_j = P(Y=y_j),
	\qqand
	M_{ij} = P(Y=y_j \mid X=x_i).
$$
Then equation \ref{eqn:cond_prob} is
\begin{align}
	P'_j = \sum_{i=1}^n P_i M_{ij}.
\end{align}
or, more compactly,
\begin{equation}
\label{eqn:P_prime_P_M}
	\vecP' = \vecP M.
\end{equation}

This matrix $M$ of transition probabilities is called a stochastic
matrix, which means the following.

\begin{defn}
A \bfidx{stochastic matrix} has all entries between 0 and 1, and all row
sums equal to 1.
\end{defn}

The reason we choose the row-sum condition is that, given $X=x_i$, it is
certain that $Y$ has to go to some $y_j$:
\begin{align*}
	\sum_{j=1}^n M_{ij}
	&=
	\sum_{j=1}^n P(Y=y_j \mid X=x_i) = 1.
\end{align*}

%% ----------------------------------------------------------------
\subsection{Markov chains}
\label{sec:defn_markov_chain}

Remember that a (discrete-time) stochastic process is just a sequence of random
variables on the same probability space $\Omega$.  Those random variables can
be related or unrelated.  We are particularly interested in the case when the
$(k+1)$st random variable depends only on the previous one.

\begin{defn}
A \bfidx{Markov chain} is a stochastic process $X^{(k)}$ such that
$$
	P(X_{n+1}=x \mid X_n=x_n, \ldots, X_1=x_1)
	=
	P(X_{n+1}=x \mid X_n=x_n).
$$
\end{defn}

The intuition is that a Markov process has (at most) only one level of memory.
The die-tipping example of section \ref{sec:die_tip} is a Markov process.  If
the up-face at the 19th step is a 2, then the up-face at the 20th step can be
1, 3, 4, or 6.  But it doesn't matter at all \emph{how} the 19th face got to be
a 2.

(I said ``at most'' since a sequence of independent random variables is also a
Markov process:  for example, replace the die-tipping experiment with a
die-rolling experiment.)

(An example of a non-Markov process is a modified die-tipping experiment
wherein I tip to a new face that avoids the previous one.  E.g.  if the 10th
face is a 3 and the 9th face was a 2, I would choose the 11th face with equal
likelihood from 1, 2, and 6.)

To apply conditional probability to Markov processes, simply replace $X$ and
$Y$ of the previous section with sequence of $W^{(k)}$ which are identity
functions from $\Omega\to\Omega$.  The $W^{(k)}$ is simply the random variable
which is the state (of the coin flips, die-tips, etc.) at the $k$th time step.
Also replace the $\vecP$'s of the previous section with a sequence of
$\vecP^{(k)}$'s.  Then we have probablity measures
$$
	P^{(k+1)}(\{\omega^{(j)}\})
	\qqand
	P^{(k)}(\{\omega^{(i)}\}).
$$
on the sample space $\Omega$, we have \plainidx{transition probabilities}
$$
	M^{(k,k+1)}_{ij} = P(W_{k+1} = \omega^{(j)} \mid W^{(k)} = \omega^{(i)}),
$$
and the PMFs at each successive time step are related by matrix multiplication:
just as in equation \ref{eqn:P_prime_P_M}:
$$
	P^{(k+1)}_j = \sum_i P^{(k)}_i M^{(k,k+1)}_{ij}
$$
i.e.
\begin{equation}
	\vecP^{(k+1)} = \vecP^{(k)} M^{(k,k+1)}.
\end{equation}
If additionally (as is the case in this paper) the transition probabilities are
the same for each $k$, then we have a \bfidx{time-homogeneous Markov chain}.
Then the successive PMFs are related by
\begin{equation}
\label{eqn:PPM}
	\vecP^{(k+1)} = \vecP^{(k)} M
\end{equation}
and so
\begin{equation}
\label{eqn:PM0}
	\vecP^{k} = \vecP^{(0)} M^k.
\end{equation}
Thus, it will be interesting below [xxx] to see what happens to powers
$M^k$ of the Markov matrix $M$.

For our purposes, think of a Markov chain as a fixed matrix $M$ along with with
the sequence $\vecP^{(k)}$ of PMFs, remembering that there's only the choice of
$\vecP^{(0)}$ after which point subsequent PMFs are obtained by matrix
multiplication.

When implementing such things in software (the electronics don't have the
error-correcting intuition your brain does), it's especially important to get
the $i$'s and $j$'s straight:
\begin{itemize}
\item $M_{ij} = P(W_{k+1} = \omega^{(j)} \mid W^{(k)} = \omega^{(i)})$.
\item $M_{ij} = P(\omega^{(i)} \mapsto \omega^{(j)})$.
\item When looking for the transition probability \emph{from} state $i$
	\emph{to} state $j$, look in row $i$ and column $j$ of $M$.
\end{itemize}

%% ----------------------------------------------------------------
\subsection{Die-tipping revisited}
\label{sec:die_tip_2}

Now that we know about the evolution of PMFs using Markov matrices, it's much
easier to compute the PMFs for the die-tipping experiments, as promised in
section \ref{sec:die_tip}.  In both experiments, the transition probabilities
are
$$
	M =
	\pcmatrix{
		0 & \frac{1}{4} & \frac{1}{4} & \frac{1}{4} & \frac{1}{4} & 0 \\
		\frac{1}{4} & 0 & \frac{1}{4} & \frac{1}{4} & 0 & \frac{1}{4} \\
		\frac{1}{4} & \frac{1}{4} & 0 & 0 & \frac{1}{4} & \frac{1}{4} \\
		\frac{1}{4} & \frac{1}{4} & 0 & 0 & \frac{1}{4} & \frac{1}{4} \\
		\frac{1}{4} & 0 & \frac{1}{4} & \frac{1}{4} & 0 & \frac{1}{4} \\
		0 & \frac{1}{4} & \frac{1}{4} & \frac{1}{4} & \frac{1}{4} & 0 \\
	}
$$
E.g. if the 1-face is up, afte the die-tip the up-face can't be 1 or 6 but it
can be any of the other four with equal likelihood.

In experiment A, the initial PMF was
$$
	\vecP^{(0)} = (1,0,0,0,0,0).
$$
Then we can just turn the crank and compute
$$
\begin{tabular}{|c|cccccc|}
\hline $k$ & \multicolumn{6}{|c|}{$\vecP^{(k)}$} \\
\hline  0 & 1.0000000 & 0.0000000 & 0.0000000 & 0.0000000 & 0.0000000 & 0.0000000 \\
\hline  1 & 0.0000000 & 0.2500000 & 0.2500000 & 0.2500000 & 0.2500000 & 0.0000000 \\
\hline  2 & 0.2500000 & 0.1250000 & 0.1250000 & 0.1250000 & 0.1250000 & 0.2500000 \\
\hline  3 & 0.1250000 & 0.1875000 & 0.1875000 & 0.1875000 & 0.1875000 & 0.1250000 \\
\hline  4 & 0.1875000 & 0.1562500 & 0.1562500 & 0.1562500 & 0.1562500 & 0.1875000 \\
\hline  5 & 0.1562500 & 0.1718750 & 0.1718750 & 0.1718750 & 0.1718750 & 0.1562500 \\
\hline  6 & 0.1718750 & 0.1640625 & 0.1640625 & 0.1640625 & 0.1640625 & 0.1718750 \\
\hline  7 & 0.1640625 & 0.1679688 & 0.1679688 & 0.1679688 & 0.1679688 & 0.1640625 \\
\hline  $\vdots$ & $\vdots$ & $\vdots$ & $\vdots$ & $\vdots$ & $\vdots$ & $\vdots$ \\
\hline 19 & 0.1666660 & 0.1666670 & 0.1666670 & 0.1666670 & 0.1666670 & 0.1666660 \\
\hline 20 & 0.1666670 & 0.1666665 & 0.1666665 & 0.1666665 & 0.1666665 & 0.1666670 \\
\hline
\end{tabular}
$$

xxx comment, and plots.

Likewise, for experiment B,

\begin{tabular}{|c|cccccc|}
\hline $k$ & \multicolumn{6}{|c|}{$\vecP^{(k)}$} \\
\hline  0 & 0.1666667 & 0.1666667 & 0.1666667 & 0.1666667 & 0.1666667 & 0.1666667 \\
\hline  1 & 0.1666667 & 0.1666667 & 0.1666667 & 0.1666667 & 0.1666667 & 0.1666667 \\
\hline  2 & 0.1666667 & 0.1666667 & 0.1666667 & 0.1666667 & 0.1666667 & 0.1666667 \\
\hline  $\vdots$ & $\vdots$ & $\vdots$ & $\vdots$ & $\vdots$ & $\vdots$ & $\vdots$ \\
\hline
\end{tabular}

xxx compute higher powers of $M$.

%% ----------------------------------------------------------------
\subsection{Markov-chain properties}

These are discussed thoroughly in chapter 4 of \cite{bib:kerldis}:
\begin{itemize}
\item \emph{irreducibility}
\item \emph{aperiodicity}
\item \emph{ergodicity}
\item \emph{reversibility}
\item \emph{stationarity}
\end{itemize}

Numerical examples from handwritten notes.

xxx plot them with and without transposing --- informative!

xxx note that this is ideal for detecting thermalization, but for large $N$,
writing down the entire PMF is precisely what we cannot do.  Hence the art of
thermalization detection.

%Given $M$, find stationary $P$:

%\begin{itemize}
%
%\item Irr and aper implies unique existence; common rows of high powers of $M$,
%with examples.
%
%\item For that, do some tables and some pygr plots.  Define thermalization
%time.
%
%\item Really not important for this paper, though.  More important:  given $P$,
%find an $M$ which has $P$ as its stationary distribution.  Xref to the section.
%
%\end{itemize}

%% ================================================================
\newpage
\section{Markov chain Monte Carlo}
\label{sec:MCMC}

%% ----------------------------------------------------------------
\subsection{Metropolis choice of Markov matrix}

Given $P$, find $M$.  Various possibilities; ``good'' convergence rate.

Def the particular $M$.

Prove it satisfies DB by construction.

Prove irr.

Prove aper.  Not true for uniform $P_0$!  Need something else there \ldots.

Note high correlation of $X^{(k)}$'s.

%% ----------------------------------------------------------------
\subsection{Sampling theory}

Thm:  Sample mean over chain converges to true mean, when $M$ is finite,
irr, and aper.

Prove this.

$c/\sqrt{n}$ convergence rate.

% [Huang section 18.9]
% Type up handwritten notes:
%
% (1) The equilibrium distribution is an eigenstate of the transition matrix.
%
% (2) The distance between subsequent PMFs is non-increasing.

%% ================================================================
\newpage
\section{Generalizations}

xxx to do.  Citations.

%% ----------------------------------------------------------------
\subsection{2D Ising}

%% ----------------------------------------------------------------
\subsection{Non-example:  percolation}

%% ----------------------------------------------------------------
\subsection{Numerical integration}

xxx continuous here \ldots.

xxx emph large number of variables; Riemann mesh impossible.

%% ================================================================
% To do:

% * re-index

% * write-up on deterministic behavior for constant-probability cases with rmk=0.
%   Affects initial state as well as rmk.

% * write-up on other random variables:  E, M, corr0j.

% * Note what the Markov matrix looks like in the independent case.

% * Display Markov matrices (with mmv option) in the independent case.

% * Display Markov matrices (with mmv option) in the 1D-Ising case.

% * sackmat assignment (copy-ctor) operator?  test out need for it.

% * plot rows of P M^k (n=3) to show thermalization

% * theoretically figure out the wide sweep ... then play with it.

% * Tight C code for larger-n sims.

% *! correlations for equal energies!!

% * therm <--> "mixing"

% * Huang:  section 18.9

% MC
%   MCMC
%     MH
%     Metro: symmetric proposal density (D.B.)
%     Hastings: allow asymmetric proposal density
%       Gibbs (Gibbs sampler); Stuart and Donald Geman

% * reversible ... how to relate to stationary?
% * symmetric mkv mx -- ?

% * Stochastic processes:  discrete-time (input) and/or discrete-valued
%   (output -- PDF/PMF).

% * find a deeper-than-Markov example.

% $\lim_{n \to \infty} P(X_n = k) = \pi(k)$.
