% This is a simple LaTex sample document that gives a submission format
%   for IEEE PAMI-TC conference submissions.  Use at your own risk.
% Make two column format for LaTex 2e.\
\documentclass[11pt,twocolumn]{article} %,twocolumn
\usepackage{times,amsmath,amsfonts}

% Use following instead for LaTex 2.09 (may need some other mods as well).
%\documentstyle[times,twocolumn]{article}
\usepackage[dvips]{graphicx,graphics}
% Set dimensions of columns, gap between columns, and paragraph indent
\setlength{\textheight}{11in} \setlength{\textwidth}{7.5in}
%\setlength{\columnsep}{0.3125in} \setlength{\topmargin}{0in}
\setlength{\headheight}{-0.5in} \setlength{\headsep}{-1in}
\setlength{\parindent}{1pc}
\setlength{\oddsidemargin}{-0.6in}  % Centers text.
\setlength{\evensidemargin}{-0.5in}

% Add the period after section numbers.  Adjust spacing.
\newcommand{\Section}[1]{\vspace{-8pt}\section{\hskip -1em.~~#1}\vspace{-3pt}}
\newcommand{\SubSection}[1]{\vspace{-3pt}\subsection{\hskip -1em.~~#1}
        \vspace{-3pt}}
\newcommand{\bqn}{\begin{eqnarray}}
\newcommand{\eqn}{\end{eqnarray}}
\newcommand {\diff}[1] {\frac{\partial}{\partial #1}}
\newcommand{\jacob}[3]{\frac{\partial^2 #3}{\partial #1 \partial #2}}
\newcommand{\der}[2]{\frac{\partial #2}{\partial #1}}
\begin{document}

% Make title bold and 14 pt font (Latex default is non-bold, 16pt)
\title{Stat 471: Lecture 21\\
MCMC I.}
% For single author (just remove % characters)
\author{Moo K. Chung\\
mchung@stat.wisc.edu}
% For two authors (default example)
\maketitle \thispagestyle{empty}
\begin{figure}
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.4]{lecture21-1.eps}
\end{figure}
\begin{enumerate} 
\item A {\em Markov chain Monte-Carlo} (MCMC) method for the simulation of a distribution $\pi$ is any method producing a Markov chain $X_i$ whose invariant distribution is $\pi$. 
\item {\em Metropolis-Hastings Algorithm} (M-H). For given distribution $\pi$ and conditional density $q(y|x)$ called {\em instrumental} or {\em proposal} distribution, define
$$\alpha(x,y) = \min \Big(\frac{\pi(y)}{\pi(x)}\frac{q(x|y)}{q(y|x)}, 1\Big).$
Then we generate Markov chain $X_i$ by\\
1. Initialize chain $X_0$ and $i=0$.\\
2. Generate $Y \sim q(y|X_i)$.\\
3. Generate $U \sim Unif(0,1)$.\\
4. If $U \leq \alpha(X_i,Y)$, $X_{i+1}=Y$\\
    else $X_{i+1} = X_i$.\\
5. Set $i = i+1$ and goto step 2.\\
Note that $0\leq \alpha(X_i,Y) \leq 1$ so the step 4 is equivalent to take $X_{i+1} = Y$ with probability $\alpha(X_i,Y)$ and $X_{i+1}=X_i$ with probability $1 - \alpha(X_i,Y)$. $q$ should be easy to simulate from to be able to use M-H algorithm. 

\item The transition kernel of the resulting chain $X_i$ is $p_{xy} = q(y|x)\alpha(x,y)$ if $y \neq x$ and $p_{xx} = 1 - \sum_{u \neq x} p_{xu}= 1- \sum_{u \neq x} q(u|x)\alpha(x,u)$. In this case it is trivial to show that 
$$\pi(x)p_{xy}=\pi(y)p_{yx} \; \mbox{(HW 4)}.$$
This is called the {\em detailed balance condition} of a Markov chain. Then it can be shown that $\pi$ is the invariant distribution and the chain is reversible, i.e. $p_{xy}=p_{yx}$.\\
 {\em Proof.} $\sum_{y \in S} \pi(y)p_{yx} = \sum_{y \in S} \pi(x)p_{xy} = \pi(x).$ From Lecture 20, this is equivalent to the invariant condition in matrix form $\pi'=\pi'{\bf P}$.

\item The stationary of $\pi$ is established for any conditional distribution $q$. This shows the universality of M-H algorithm. Based on M-H algorithm and the ergodic theorem, then we can estimate $\mathbb{E}_{\pi} g \approx \frac{1}{n-m}\sum_{i=m+1}^n g(X_i)$, where the initial chain $X_0$ and the conditional distribution $q$ do not really matter.

\item For large $m$, $X_m, X_{m+1}, \cdots$ produces a dependent sample distributed from $\pi$ (see Lecture 20 for detail). Let's generate 10000 Cauchy numbers from normal distribution using MCMC. Note that they will not be independent. Let $\pi(x) = \frac{1}{1+x^2}$ (no need to specify the normalizing constant) and $q(y|x) \sim N(x,1)$

\begin{verbatim}
n=10000;
x(1)=0;
str='exp(-(x-mu)^2/2)';
snrpdf=inline(str,'x','mu');
cauchypdf=inline('1/(1+x^2)');
for i=2:n
  y=x(i-1) + snrnd(1);
  alpha=min([1,snrpdf(y)/snrpdf(x(i-1))
  *snrpdf((x(i-1),y)/snrpdf(y,x(i-1))]);
  u=rand;
  if u <= alpha
    x(i)=y;
  else
    x(i)=x(i-1);
  end 
end  
\end{verbatim}
\begin{figure}
\centering
\renewcommand{\baselinestretch}{1}
\includegraphics[scale=0.4]{lecture21-2.eps}
\end{figure}
\end{enumerate} 
\end{document}

n=10000;
x(1)=0;
str='exp(-(x-mu)^2/2)';
snrpdf=inline(str,'x','mu');
cauchypdf=inline('1/(1+x^2)');
for i=2:n
  y=x(i-1) + randn;
  alpha=min([1,cauchypdf(y)/cauchypdf(x(i-1))*snrpdf(x(i-1),y)/snrpdf(y,x(i-1))]);
  u=rand;
  if u <= alpha
    x(i)=y;
  else
    x(i)=x(i-1);
  end 
end