% Borrowed from MIT 6.987 Advanced Data Structures course (Prof Demaine, 2003)

\documentclass[11pt]{article}
\usepackage[utf8]{inputenc} 
\usepackage{amsmath, amssymb, amsthm, bm}
\usepackage[usenames,dvipsnames]{xcolor}
\usepackage[colorlinks,citecolor=blue,urlcolor=blue,linkcolor=blue]{hyperref}
\usepackage{hypernat}

\usepackage{tikz}
\usepackage{tkz-graph}

\usepackage{listings}             % Include the listings-package
\usepackage{algorithm2e}

\newcommand{\handout}[5]{
  \noindent
  \begin{center}
  \framebox{
    \vbox{
      \hbox to 5.78in { {\bf STA4516: Topics in 
        Probabilistic Programming } \hfill #2 }
      \vspace{4mm}
      \hbox to 5.78in { {\Large \hfill #5  \hfill} }
      \vspace{2mm}
      \hbox to 5.78in { {\em #3 \hfill #4} }
    }
  }
  \end{center}
  \vspace*{4mm}
}

\newcommand{\lecture}[4]{\handout{#1}{#2}{#3}{Scribe: #4}{Lecture #1}}

\newtheorem{theorem}{Theorem}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{observation}[theorem]{Observation}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{claim}[theorem]{Claim}
\newtheorem{fact}[theorem]{Fact}
\newtheorem{assumption}[theorem]{Assumption}
\newtheorem{remark}[theorem]{Remark}

% 1-inch margins, from fullpage.sty by H.Partl, Version 2, Dec. 15, 1988.
\topmargin 0pt
\advance \topmargin by -\headheight
\advance \topmargin by -\headsep
\textheight 8.9in
\oddsidemargin 0pt
\evensidemargin \oddsidemargin
\marginparwidth 0.5in
\textwidth 6.5in

\parindent 0in
\parskip 1.5ex
%\renewcommand{\baselinestretch}{1.25}

\begin{document}

\lecture{2 --- November 3, 2015}{Fall 2015}{Prof.\ Daniel M. Roy}
{Mufan (Bill) Li}

\section{Overview}

In the last lecture we saw an introduction and motivations
for studying probabilistic programs.

This lecture partly taught using presentation slides and 
partly on chalk board.
In the presentation, we briefly discussed 
\textbf{exchangeable} structures, 
\textbf{de Finetti's Theorem} and the application to achieve 
\textbf{conditional independence}, 
and the \textbf{computability} properties of these structures.
On the chalk board, we saw a lightweight implementation of 
\textbf{probabilistic programs} 
using a \textbf{Metropolis-Hasting} type algorithm.

\section{Presentation}

Here I add several notes supplementary to the 
presentation,
please see the original slides for full details.

\subsection{de Finetti's Theorem}

An exchangeable sequence $(Y_1, Y_2, \ldots, Y_N)$
has a fully connected graphical model, 
which has $\mathcal{O}(N^2)$ number of connections.
By applying de Finetti's Theorem, 
we are able to reduce the number of connections to $N$,
as well as allow \textbf{parallel computation} of random numbers
since we have conditional independence.

An important remark is the transformation $f(\theta,U)$,
where conditioned on $\theta$ gives independence and 
$U$ is a uniform random variable,
is \textbf{not unique}.
Here any measure-preserving 
transformations of $\theta$ and $U$ 
will give the same results,
specifically
%
\begin{equation*}
  \psi(U) \,{\buildrel d \over =}\ U \Rightarrow
  \exists f^\psi(\theta, \psi(U)):
  \big(Y_1, \dots, Y_N\big) \,{\buildrel d \over =}\
  \big(f^\psi(\theta, \psi(U_1)), \ldots, f^\psi(\theta, \psi(U))\big)
\end{equation*}

\section{Chalk Board}

In this section, we discuss a lightweight implementation
of probabilistic programs from \cite{WSG11}.

\subsection{Introduction}

\begin{definition}
  An \textbf{unconditional probabilistic program}
  is a parameterless function $f$ with an arbitrary mix of 
  stochastic and deterministic elements.
\end{definition}

\textbf{Informally} $f$ is a deterministic program
until it calls an external stochastic coroutine
to generate randomness.
For example, MATLAB calls \textit{rand} to generate 
a uniform random variable.
These coroutines are called 
\textbf{elementary random primitives} (ERPs).

\begin{definition}
  The families of distributions
  $\big(P_t(\cdot, \theta)\big)_t$ 
  are \textbf{dominated} by some $\sigma$-finite
  measure $\mu_t$ if
  %
  \begin{equation*}
    \exists \text{ densities } p_t(\cdot, \theta) :
    P_t(A|\theta) = \int_A p_t(x|\theta) \mu_t(dx)
  \end{equation*}
\end{definition}

For example, the family of Gaussian distributions
has densities $p_t(\cdot|\theta) = \mathcal{N}(\overline{x}, \sigma^2)$,
where $t$ determines the Gaussian type, 
and $\theta$ determines the parameters (mean and variance)
$(\overline{x}, \sigma^2)$.
Here the Gaussian family is dominated by the Lebesgue
measure i.e $\mu_t(dx) = m(dx)$.

\begin{definition}
  $\mathcal{T}$ is a \textbf{set of ERP types}, 
  where for each $t \in \mathcal{T}$, 
  we have a dominated family of distributions
  $\big(P_t(\cdot, \theta)\big)_t$, with densities
  $p_t(x | \theta)$ where $\theta \in \Theta_t$.
\end{definition}

Note $\mu_t$ can depend on the ERP type $t$ 
but not parameters $\theta$. 
It is also important to note the requirement that these
families of distributions to be dominated is rarely mentioned,
but implicit in the way that these densities are used.

Here we introduce a simple example in MATLAB.

\lstset{language=Matlab}
\begin{lstlisting}[frame=single]  % Start your code-block

for i = 1:1000
  if (rand > 0.5)
    x(i) = randn;
  else
    x(i) = gammarnd;
  end
end
\end{lstlisting}
%
where \textit{rand} generates a uniform random variable,
\textit{randn} corresponds to a standard normal,
and \textit{gammarnd} is gamma.

Depending on the return values of ERPs, 
we may have different paths of execution.
Here we introduce the \textbf{execution trace},
a minimal log of execution containing 
the requests for ERP and return values,
which is sufficient to reconstruct the path of execution.

A possible trace may look like this:

\begin{tabular}{ l | l | l | l }
  Index & $t$ & $\theta$ & Return \\
  \hline
  1 & rand & $[a,b]$ & 0.7 \\
  2 & randn & $(\overline{x},\sigma^2)$ & -2.7 \\
  3 & rand & $[a,b]$ & 0.2 \\
  4 & gammarnd & $(k,\phi)$ & 1.9 \\
  $\vdots$ & $\vdots$ & $\vdots$ & $\vdots$\\
  2000 & randn & $(\overline{x},\sigma^2)$ & 1.5 \\
\end{tabular}

Here a total of 2000 random choices are made.
It is also possible to for it to be random.
For example:

\begin{lstlisting}[frame=single]  % Start your code-block

accept = false;
n = 0;
list = [];
while not accept
  n = n + 1;
  list = [list n]; % concatenate n to the end of list
  accept = rand < 0.5
end
\end{lstlisting}

The number of ERPs encountered here is geometrically 
distributed with mean 2.


\subsection{Statistical Inference as Conditional Simulation}

Consider the following program

\begin{lstlisting}[frame=single]  % Start your code-block

X = gammarnd;
for j = 1:1000
  Y(j) = poissonrnd(X);
end
\end{lstlisting}

This program generates a sequence of Poisson random variables
using a gamma prior, \\
i.e. $X \sim \text{ Gamma}, \; 
Y_j \sim \text{ Poisson}(X^{-1}), Y = (Y_1,\ldots,Y_{1000}) $.
If $Y_j$ model the number of emails that employee $j$ 
receives in one day,
how would we estimate the (common) rate at which 
email arrives?

Suppose we had data $y = (y_1,\ldots,y_{1000})$, 
we could simulate the program many times until we see
data $Y = y$.
The value $X$ used for that particular simulation
is then a perfect sample from the posterior
$P(X|Y = y)$. This is also known as \textbf{rejection sampling},
\textbf{accept-reject}, or \textbf{guess and check}.

By \textbf{conditional simulation}, 
we mean a simulation of the program sampled
from the distribution over simulations satisfying
some predicate.
In this case the predicate is ``accepts'' a simulation
if $Y = y$, and rejects it otherwise.

\subsubsection{Computing Posterior Distributions}

In general sampling and taking expectations
from a conditional distribution 
$P(X|Y=y)$ can be difficult even if we can calculate the 
joint density $p(x,y)$.
The reason is conditional densities has the form 
%
\begin{equation*}
  p(x|y) = \frac{p(x,y)}{p(y)} = \frac{p(x,y)}{\int p(x,y) 
    \tau(dx \times Y)}
\end{equation*}
%
and the integral on the bottom is generally hard to compute,
even when given $p$ and $\tau$.
Even if the conditional density is known,
to sample efficiently we also need the following
%
\begin{equation*}
  P(x \in A|Y = y) = \int_A p(x|y) \tau(dx \times Y)
\end{equation*}
%
to be computed efficiently as well.
In order to address this issue, we consider the
Metropolis-Hasting algorithm.

\subsubsection{Metropolis-Hasting (MH) Algorithm}

MH is a way to sample from the posterior $P(X|Y=y)$ 
or compute conditional
expectations $\mathbb{E}[f(X)|Y=y]$.
MH is one of many 
\textbf{Markov chain Monte Carlo} (MCMC) algorithms.
For an excellent survey of MCMC algorithms see \cite{RR04}.

We first fix a ``proposal distribution'',
or more formally a 
\textbf{transition probability kernel} \\
%
\begin{equation*}
  Q : X \to M_1(X)
\end{equation*}
%
on the space $X$, where $M_1(X)$ is the space of 
distributions on $X$.
In other words,
$Q$ is a family of distributions parameterized by $x \in X$.

We are going to assume that $Q$ is dominated, i.e.
there exist a $\sigma$-finite measure $\nu$ 
and density $q$ such that 
%
\begin{equation*}
  Q(x,A) = Q(x)(A) = \int_A q(x,y) \nu(dy)
\end{equation*}
%
For example, a choice of $Q$ for random walk 
Metropolis-Hasting (RWMH) is the following
%
\begin{equation*}
  \epsilon \sim \mathcal{N}(0,1), \;
  Q(x) = \mathcal{L}(x+\epsilon) \,{\buildrel d \over =}\
  \mathcal{N}(x,1)
\end{equation*}
%
Here $Q(x)$ is a measure, $q(x,y)$ is Gaussian,
and $\nu(dy)$ is the Lebesgue measure.
This leads to a small perturbation of current value of $x$.
Note the choice of $\nu$ is important if the space of 
interest can be a null set w.r.t. $\nu$.
This is possible if the space of interest is lies on  
a sub-manifold of measure zero, such as a curve in $\mathbb{R}^2$.

We then construct the following Markov chain

\begin{algorithm}[H]
 % \KwData{this text}
 % \KwResult{how to write algorithm with \LaTeX2e }
 Initialize $X_0, i=1$\;
 \Repeat{}{
  Sample $X_i' \sim Q(X_{i-1})$\;
  \KwSty{with probability} $a(X_{i-1},X_i)$ \\
  \Indp
  Set $X_i = X_i'$ (accept)\;
  \Indm
  \KwSty{otherwise}\\
  \Indp
  Set $X_i = X_{i-1}$ (reject)\;
  \Indm
  $i = i+1$\;
 }
 % \caption{How to write algorithms}
\end{algorithm}
%
where we define the function $a$ as 
%
\begin{equation*}
  a(x,x') = \min\left\{1, \frac{p(x'|y)q(x',x)}{p(x|y)q(x',x)}
    \right\}
\end{equation*}

The choice of the $a(x,x')$ results in the following property:
let $T(x,\cdot) = P(X_1|X_0 = x)$ be the transition
probability kernel for the chain $(X_0, X_1, \ldots)$, then 
%
\begin{equation}\label{eq:detailed_balance}
  P(X \in dx | Y=y) T(x, dx') = P(X \in dx'|Y=y) T(x',dx)
\end{equation}
%
Therefore the Markov chain is \textbf{reversible}.
In other words, the chain $(X_0, X_1, \ldots)$
can be extended to $(\ldots, X_{-2}, X_{-1}, X_0)$
and the backward processes are equal in distribution.

Condition \eqref{eq:detailed_balance} is called 
\textbf{detailed balance}.
When \eqref{eq:detailed_balance}
holds, it implies that $P(X|Y=y)$ is a stationary distribution:
i.e. if we sample $X_0 \sim P(X|Y=y)$, 
then every $X_j \sim P(X|Y=y)$.
Furthermore, if the Markov chain is \textbf{ergodic}, 
then for an arbitrary starting point $X_0$, we have
%
\begin{equation*}
  X_j \xrightarrow{d} X \sim P(X|Y=y)
\end{equation*}
%
That is, we can generate a sample from the posterior 
$P(X|Y=y)$ from running the Markov chain from 
any starting point $X_0$.

\begin{fact}
  Only requirement of $Q$ is ergodicity. Here the choice
  of $a$ will bias $Q$ to achieve detailed balance
  \eqref{eq:detailed_balance}.
\end{fact}

\begin{observation}
  The MH algorithm allowed us to avoid computing
  the integral by taking a ratio instead
  %
    \begin{equation*}
      \frac{P(x'|y)}{P(x|y)} = \frac{P(x',y)}{P(x,y)}
    \end{equation*}
  %
  This is a result of ergodicity, 
  where we exchange the integral over space to 
  integral over time.
\end{observation}

\subsubsection{MH in Probabilistic Programming}

In order to apply MH, we need to answer the 
following questions:
%
\begin{enumerate}
  \item What is the space $X$? \\
    Space of trace.
  \item What is the transition probability kernel $Q$?\\
    Transition probability from one trace path to another.
  \item How to compute $a$? \\
    TBA next lecture.
\end{enumerate}
%
\begin{thebibliography}{9}

\bibitem{WSG11}
  David Wingate, Andreas Stuhlmüller, and Noah D. Goodman,
  \emph{Lightweight Implementations of Probabilistic 
  Programming Languages Via Transformational Compilation},
  AISTATS, 2011.
\bibitem{RR04}
  Gareth O. Roberts, Jeffrey S. Rosenthal,
  \emph{General state space Markov chains and MCMC algorithms}.
  Probability Surveys, 2004.
\end{thebibliography}
\bibliographystyle{alpha}

\end{document}
