\documentclass[11pt, fleqn]{article}
%\VignetteIndexEntry{GiRaF-introduction}
%\VignetteEngine{knitr::knitr}
%\VignetteKeywords{GiRaF, NC.mrf, exact.mrf, sampler.mrf, mrf, grf, potts, Gibbs}
%\VignetteDepends{GiRaF}

%-------------------------
% Gestion des langues
%-------------------------
\usepackage[utf8]{inputenc}
\usepackage[english]{babel}

%-------------------------
% Gestion des fontes
%-------------------------
\usepackage{fourier} 
%\usepackage{pifont}

%-------------------------
% Gestion des packages
%-------------------------

% For bibliography
\usepackage{natbib}
\bibpunct{(}{)}{,}{a}{,}{,}

% For maths
\usepackage{amsmath}
\makeatletter
\def\resetMathstrut@{%
  \setbox\z@\hbox{%
    \mathchardef\@tempa\mathcode`\(\relax
      \def\@tempb##1"##2##3{\the\textfont"##3\char"}%
      \expandafter\@tempb\meaning\@tempa \relax
  }%
  \ht\Mathstrutbox@1.2\ht\z@ \dp\Mathstrutbox@1.2\dp\z@
}
\makeatother
\usepackage{amssymb}
\usepackage{amsfonts}
\usepackage{amsthm}

\usepackage{cancel}
\usepackage{mathrsfs}
\usepackage{dsfont}

% For pictures
%\graphicspath{{images/}}
\usepackage{graphicx}
\usepackage[format=hang, font=small, labelfont=bf]{caption}
\usepackage{tikz} 


% For diagrams
%\usepackage[all]{xy}

% For table
\usepackage{array}
\usepackage{colortbl}
\usepackage{multicol}
\usepackage{color}
\definecolor{mygreen}{rgb}{0.16,0.5,0.32}
\definecolor{gris50}{gray}{0.50}
\definecolor{gris95}{gray}{0.95}
\definecolor{myblue}{rgb}{0.22,0.45,0.70}


% For algorithm
\usepackage[algo2e,ruled]{algorithm2e}

% For hyperref
\usepackage[a4paper, citecolor = blue,
            linkcolor = black,
            colorlinks = true, 
            urlcolor=blue]{hyperref}


% Page layout
\usepackage{enumitem}
\usepackage{chngpage}
\usepackage{pdfpages}

\usepackage{geometry}
\setlength{\textwidth}{146.8mm} % = 210mm - 37mm - 26.2mm
\setlength{\oddsidemargin}{11.6mm} % 37mm - 1in (from hoffset)
\setlength{\evensidemargin}{0.8mm} % = 26.2mm - 1in (from hoffset)
\setlength{\topmargin}{-2.2mm} % = 0mm -1in + 23.2mm 
\setlength{\textheight}{221.9mm} % = 297mm -29.5mm -31.6mm - 14mm (12 to accomodate footline with pagenumber)
\setlength{\headheight}{14pt}%{14.49998pt}

\setlength{\parindent}{0pt}
\usepackage{setspace} % increase interline spacing slightly
\setstretch{1.1}
\usepackage{parskip}

\usepackage{fancyhdr}
\pagestyle{fancy}

\renewcommand{\headrulewidth}{1pt}
\fancyhead[L]{J. Stoehr, P. Pudlo and N. Friel}
\fancyhead[R]{GiRaF: a toolbox for Gibbs Random Field Analysis}

\usepackage{listings}
\lstset{language=[LaTeX]Tex,tabsize=4, basicstyle=\scriptsize\ttfamily, showstringspaces=false, numbers=left, numberstyle=\tiny, numbersep=10pt, breaklines=true, breakautoindent=true, breakindent=10pt}


\usepackage{authblk}

\title{\textbf{GiRaF: a toolbox for Gibbs Random Fields analysis}}
%\title{\textcolor{myblue}{\textbf{GiRaF: a toolbox for Gibbs Random Fields analysis}}}
\author[1]{Julien Stoehr\thanks{julien.stoehr@ucd.ie}}
\author[2]{Pierre Pudlo}
\author[1]{Nial Friel}
\affil[1]{University College Dublin}
\affil[2]{Aix-Marseille Université}
\date{\today}


%-------------------------
% Commandes perso
%-------------------------

\newtheoremstyle{Rstyle}
  {5pt} % Space above
  {\topsep} % Space below
  {} % Body font
  {} % Indent amount
  {\bfseries} % Theorem head font
  {:} % Punctuation after theorem head
  {.5em} % Space after theorem head
  {} % Theorem head spec (can be left empty, meaning `normal')
\newtheoremstyle{Rexample}
  {20pt} % Space above
  {\topsep} % Space below
  {} % Body font
  {} % Indent amount
  {\bfseries} % Theorem head font
  {:} % Punctuation after theorem head
  {.5em} % Space after theorem head
  {} % Theorem head spec (can be left empty, meaning `normal')
\theoremstyle{Rstyle}
\newtheorem*{inputs}{Inputs}
\newtheorem*{usage}{Usage}
\theoremstyle{Rexample}
\newtheorem*{example}{Example}

\newcommand{\suchthat}{\;\ifnum\currentgrouptype=16 \middle\fi|\;}

\newcommand{\s}{\mathscr{S}}
\newcommand{\G}{\mathscr{G}}

\newcommand{\Xv}{\mathbf{X}}
\newcommand{\x}{\mathbf{x}}
\newcommand{\X}{\mathscr{X}}
\newcommand{\Y}{\mathscr{Y}}
\newcommand{\y}{\mathbf{y}}
\newcommand{\Yv}{\mathbf{Y}}
\newcommand{\yobs}{\mathbf{y}^{\text{obs}}}
\newcommand{\xobs}{\mathbf{x}^{\text{obs}}}

\newcommand{\ivj}{i\stackrel{\G}{\sim} j}
\newcommand{\lvj}{\ell\stackrel{\G}{\sim} j}
\newcommand{\ivjf}{i\stackrel{\G_4}{\sim} j}

\newcommand{\N}{\mathscr{N}}
\newcommand{\C}{\mathscr{C}}

\newcommand{\mle}{\text{MLE}}
\newcommand{\mcle}{\text{MCLE}}
\newcommand{\obs}{\text{obs}}
\newcommand{\BIC}{\text{BIC}}
\newcommand{\PLIC}{\text{PLIC}}
\newcommand{\BLIC}{\text{BLIC}}
\newcommand{\BCLIC}{\text{BCLIC}}
\newcommand{\mfl}{\text{MF-like}}
\newcommand{\mf}{\text{MF}}
\newcommand{\gbf}{\text{GBF}}
\newcommand{\vem}{\text{VEM}}
\newcommand{\map}{\text{MAP}}
\newcommand{\mapcl}{\text{CL}}
\newcommand{\mapclcal}{\text{cal}}
\newcommand{\thetamv}{\hat{\theta}^{\text{ml}}}

%% Proba
\newcommand{\pr}{\mathbf{P}}
\newcommand{\esp}{\mathbf{E}}
\newcommand{\ind}{\mathbf{1}}
\newcommand{\var}{\textbf{Var}}
\newcommand{\cov}{\textbf{Var}}
%\newcommand{\argmin}{\operatorname{arg\,min}}
%\newcommand{\argmax}{\operatorname{arg\,max}}
\DeclareMathOperator*{\argmin}{arg\,min\,}
\DeclareMathOperator*{\argmax}{arg\,max\,}
\DeclareMathOperator*{\gO}{\mathcal{O}}
\DeclareMathOperator*{\grad}{\mathbf{\nabla}}
\DeclareMathOperator*{\hess}{\mathbf{\nabla}^2}
%\DeclareMathOperator*{\dim}{dim}
%\DeclareMathOperator*{\trace}{tr}

%% Composite
\newcommand{\A}{\mathscr{A}}
\newcommand{\B}{\mathscr{B}}
\newcommand{\fcl}{f_{\text{CL}}^{\text{cal}}}
\newcommand{\fclu}{f_{\text{CL}}}
\newcommand{\pCL}{\pi^{\text{cal}}_{\text{CL}}}
\newcommand{\cpCL}{\overline{\pi}^{~\rm cal}_{\text{CL}}}
\newcommand{\pCLu}{\pi_{\text{CL}}}
\newcommand{\CL}{\text{CL}}
\newcommand{\trace}{\text{tr}}
\newcommand{\determinant}{\text{det}}

\newcommand{\vertiii}[1]{{\left\vert\kern-0.25ex\left\vert\kern-0.25ex\left\vert #1 
    \right\vert\kern-0.25ex\right\vert\kern-0.25ex\right\vert}}
\newcommand{\KL}{\text{KL}} 

\newcommand{\nobs}{n_{\text{obs}}}

%ABC
\newcommand{\M}{\mathscr{M}}
\newcommand{\D}{\mathcal{D}}
\newcommand{\Dobs}{\mathcal{D}^{obs}}
\newcommand{\Ss}{\mathbf{S}}

\newcommand{\p}{\pi}%{\mathbb{P}}

\newcommand{\Z}{\mathbf{Z}}
\newcommand{\indxixj}{\mathds{1}\{x_i=x_j\}}
\newcommand{\indxyi}{\mathds{1}\{x_i=y_i\}}
\newcommand{\nsamp}{N}%{n_{\text{SAMP}}}
\newcommand{\nref}{n_{\text{REF}}}
\newcommand{\npost}{n_{\text{POST}}}
\newcommand{\nnei}{n_{\text{NEI}}}
\newcommand{\HPM}{\text{HPM}(\G,\phi,\beta)}
\newcommand{\HPMf}{\text{HPM}(\G_4,\phi,\beta)}
\newcommand{\HPMh}{\text{HPM}(\G_8,\phi,\beta)}
\newcommand{\Ld}{\text{L}^2}
\newcommand{\Card}{\text{Card}}
\newcommand{\dd}{\mathrm{d}}

\newcommand{\ie}{\textit{i.e.}}
\newcommand{\eg}{\textit{e.g.}}
\newcommand{\apost}{\textit{a posteriori}~}
\newcommand{\apriori}{\textit{a priori}~}
\newcommand{\giraf}{\textbf{GiRaF}~}

\hyphenation{a-ni-so-tro-pic}

<<echo = FALSE>>=
library(GiRaF)
library(knitr)
options(prompt = " ", continue='  ', width=60)
@

\begin{document}
\maketitle

%\hrule
\begin{abstract}
\giraf package offers various tools for the analysis of Gibbs (or discrete Markov) random fields. The latter form a class of intractable statistical models since, due to the Markovian dependence structure, the normalising constant of the fields cannot be computed using standard analytical or numerical methods. \giraf substantially lowers the barrier for practitioners aiming at analysing such Gibbs random fields. It contains exact methods for small lattices and several approximate methods for larger lattices that make the analysis easier for practitioners. This document provides a short description of the models which are solved as well as short examples of the major functionalities.
\end{abstract}

%\hrule
\tableofcontents

\section*{Introduction}
Gibbs (or discrete Markov) random fields have been used in many practical settings, surged by the development in the statistical community since the 1970's. Notable examples are the autologistic model \citep{besag74} and its extension the Potts model used to describe the spatial dependency of discrete random variables (\eg, shades of grey or colors) on the vertices of an undirected graph (\eg, a regular grid of pixels). Despite their popularity, running a statistical analysis is difficult due to the dependence structure that leads to a considerable computational curse: the normalising constant is intractable and cannot be computed with standard analytical or numerical methods. 

The \giraf package aims at providing a toolbox for practitioners that makes the analysis of such models easier. The primary tool is the computation of the intractable normalising constant for small rectangular lattices (\nameref{ncmrf}) based on the recursive algorithm of \cite{reeves2004}. In practice, this function allows to have ground truth against which to compare. Beside the latter function, the \giraf package contains methods that give exact sample from the likelihood for small enough rectangular lattices (\nameref{exactmrf}, \cite{friel2007}) or approximate samples from the likelihood using MCMC samplers: the Gibbs sampler or the Swendsen-Wang algorithm (\nameref{samplermrf}) for large lattices.

The \giraf package provides \textsf{R} functions based on \textsf{C++} code for the statistical analysis of Gibbs random fields. As befits an introduction to the \giraf package, a brief overview of Gibbs distribution is presented in Section \ref{sec:markov-gibbs} while pointing out the convention used throughout the package. Then we present the major functions through detailed examples, namely recursive computing in Section \ref{sec:recursion} and sampling methods in Section \ref{sec:sampler}.

\section{Gibbs Random Fields}
\label{sec:markov-gibbs}

\subsection{Gibbs-Markov equivalence}
A discrete random field $\Xv$ is a collection of random variables $X_i$ indexed by a finite set  
$\s = \{1,\dots,n\}$, whose elements are called sites, and taking values in a finite state space 
$\mathscr{X}_i:=\{0,\ldots,K-1\}$, interpreted as colors. For a given subset $A\subset\s$, $\Xv_A$ 
and $\x_A$ respectively define the random process on $A$, \ie, $\{X_i, i\in A\}$, and a realisation
of $\Xv_A$. Denotes $\s\setminus A= -A$ the complement of $A$ in $\s$.
When modeling a digital image, the sites lie on a
regular 2D-grid of pixels, and their dependency is given by an
undirected graph $\G$ which induces a topology on $\s$: by definition, sites $i$ and $j$ are adjacent or are neighbor of each other if and only if $i$ and $j$ are linked by an edge in $\G$. A random field $\Xv$ is a Markov random field with respect to $\G$, if for all configuration $\x$ and for all sites $i$
\begin{equation}
\pr\left(X_i = x_i\suchthat \Xv_{-i} = \x_{-i}\right) = \pr\left(X_i = x_i\suchthat \Xv_{\N(i)} = \x_{\N(i)}\right),
\label{eqn:markov}
\end{equation}
where $\N(i)$ denotes the set of all the adjacent sites to $i$ in $\G$.
The Hammersley-Clifford theorem states that if the distribution of a Markov random field with respect to a graph $\G$ is positive for all configuration $\x$ then it admits a Gibbs representation for the same topology (see for example \cite{grimmett1973, besag74} and for a historical perspective \cite{clifford1990}),  namely a probability measure $\pi$ on $\X=\prod_{i=1}^n\X_i$ given by
\begin{equation}
\pi\left(\x\suchthat \psi, \G\right) = \frac{1}{Z\left(\psi, \G\right)}\exp\left\lbrace- H\left(\x\suchthat\psi,\G\right)\right\rbrace,
\label{eqn:gibbs}
\end{equation}
where $\psi$ is a vector of parameters and $H$ denotes the energy function or Hamiltonian. 

The \giraf package solely focuses on Potts model, that is a pairwise Markov random field whose Hamiltonian linearly depends on the parameter $\psi=\left\lbrace\alpha,\beta\right\rbrace$, where $\alpha$ stands for the parameter on sites and $\beta$ stands for the parameter on edges. More precisely,
\begin{equation}
H\left(\x\suchthat\alpha,\beta,\G\right) = -\sum_{i=1}^n \sum_{k=0}^{K-1}\alpha_k \ind\{x_i=k\} - \sum_{\ivj} \beta_{ij}\ind\{x_i=x_j\},
\label{eqn:hamiltonian}
\end{equation}
% \[
% H\left(\x\suchthat\psi,\G\right) = -\psi^{T}\Ss(\x).
% \]
%where $\Ss(\x)=(s_1(\x),\ldots,s_d(\x))$ is a vector of sufficient statistics. 
where the above sum $\sum_{\ivj}$ ranges the set of edges of the graph $\G$. 

At that stage, the \giraf package handles homogeneous\footnote{Parameters $\alpha$ and $\beta$ are independent of the relative position of sites or edges} and potentially anisotropic\footnote{Parameter $\beta$ depends on the orientation of the edges.} Potts models defined on a rectangular $n = h\times w$ lattice, with $h\leq w$. Lattice points are ordered from top to bottom in each column and columns from left to right. The \giraf package contains two widely used adjacency structures, namely the graph $\G_4$ (first order lattice), respectively $\G_8$ (second order lattice), for which the neighborhood of a site is composed of the four, respectively eight, closest sites on a two-dimensional regular lattice, except on the boundaries of the lattice, see Figure \ref{fig:neigh}.
\begin{figure}[t]
\centering
\begin{minipage}[t]{7cm}
\centering
\input{g4.tex}\\
(a)
%\caption{Graph $\G_4$}
\end{minipage}
\begin{minipage}[t]{7cm}
\centering
\input{g8.tex}\\
(b)
%\caption{Graph $\G_8$}
\end{minipage}
\caption{\small First and second order neighborhood graphs $\G$. (a) The four closest neighbors graph $\G_4$. Neighbors of the vertex in black are represented by vertices in gray. (b) The eight closest neighbors graph $\G_8$. Neighbors of the vertex in black are represented by vertices in gray.}
\label{fig:neigh}
\end{figure}

Within the set of allowed distributions, the default model in \giraf is set to the 2-state Potts models with a first order dependency structure $\G_4$, namely
\[
\pi\left(\x\suchthat\beta,\G_4\right) = \frac{1}{Z\left(\beta,\G_4\right)}\exp\left(\sum_{\ivjf}\beta_{ij}\ind\{x_i=x_j\}\right).
\]

\subsection{GiRaF inputs}
\label{subsec:inputs}

Most of the functions in the \giraf package rely on the following inputs which are detailed throughout this Section.
\begin{center}
\begin{tabular}{p{0.20\textwidth}p{0.74\textwidth}}
\textsf{\bf h, w} & respectively the number of rows and columns of the rectangular lattice $\G$. \\
\textsf{\bf param} & \textsf{numeric} entry setting the interaction parameter $\beta$ (edges parameter). \hyperref[par:param]{See details below}.\\
\textsf{\bf ncolors} & the number $K$ of states (or colors) for the discrete random variables $\Xv_i$. By default, $\textsf{ncolors}=2$. \\
\textsf{\bf nei} & the number of neighbors. The latter must be one of $\textsf{nei}=4$ or $\textsf{nei}=8$, which respectively correspond to $\G_4$ and $\G_8$. By default, $\textsf{nei}=4$.\\
& \\
\multicolumn{2}{l}{Optional inputs} \\
\textsf{\bf pot} & \textsf{numeric} entry setting homogeneous potential $\alpha$ on singletons. By default, \textsf{pot = NULL}.\\
\textsf{\bf top, left, bottom, right, corner} & \textsf{numeric} entry setting constant borders for lattice $\G$. By default,
\textsf{top = NULL}, \textsf{left = NULL}, \textsf{bottom = NULL}, \textsf{right = NULL}, \textsf{corner = NULL}. \hyperref[par:borders]{See details below}.\\
\end{tabular}
\end{center}

\paragraph*{Interaction parameter (\textsf{param})\label{par:param}} The \giraf packages allows to set an isotropic or anisotropic interaction parameter $\beta$ via the entry \textsf{param}. In the isotropic configuration, \textsf{param} is simply a scalar while in the anisotropic case \textsf{param} is a vector $\textsf{c}(\beta_1, \beta_2)$ or $\textsf{c}(\beta_1, \beta_2, \beta_3, \beta_4)$ if $\G=\G_4$ or $\G=\G_8$, respectively.
The parameter $\beta_{ij}$ in \eqref{eqn:hamiltonian} stands for the component of \textsf{param} corresponding to the direction defined by the edge $(i,j)$ but does not depend on the actual position of sites $i$ and $j$, that is, given two edges $(i_1,j_1)$ and 
$(i_2,j_2)$ defining the same direction, $\beta_{i_1,j_1} = \beta_{i_2,j_2}$ (see Table \ref{tab:param}).

\newcolumntype{M}[1]{>{\centering\arraybackslash}m{#1}}
\begin{table}%[b!]
\caption{\small Interaction parametrisation for a homogeneous Gibbs random field for the isotropic and anisotropic cases. The table gives values of the parameter $\beta_{ij}$ corresponding to the orientation of the edge $(i,j)$.}
\centering
\begin{tabular}{M{3cm}|M{0.75cm}M{0.75cm}|M{0.75cm}M{0.75cm}|M{0.75cm}M{0.75cm}|M{0.75cm}M{0.75cm}}

Orientation of edge $(i,j)$ 
& \multicolumn{2}{M{1.5cm}}{~~~~\includegraphics [height = 1cm, width = 1cm]{clique_v.png}} 
& \multicolumn{2}{M{1.5cm}}{~~~\includegraphics [height = 1cm, width = 1cm]{clique_h.png}}
& \multicolumn{2}{M{1.5cm}}{~~\includegraphics [height = 1cm, width = 1cm]{clique_tr.png}}
&  \multicolumn{2}{M{1.5cm}}{~~\includegraphics [height = 1cm, width = 1cm]{clique_br.png}} \tabularnewline
\hline \hline
\rowcolor[gray]{0.85} Dependency graph 
& $\G_4$ & $\G_8$
& $\G_4$ & $\G_8$
& \multicolumn{2}{M{1.5cm}|}{$\G_8$}
& \multicolumn{2}{M{1.5cm}}{$\G_8$} \tabularnewline

Isotropic Gibbs 
& \multicolumn{8}{c}{$\beta$} \tabularnewline
\rowcolor[gray]{0.85} Anisotropic Gibbs
& \multicolumn{2}{c|}{$\beta_1$}
& \multicolumn{2}{c|}{$\beta_2$}
& \multicolumn{2}{c|}{$\beta_3$}   
& \multicolumn{2}{c}{$\beta_4$}
\label{tab:param}
\end{tabular}
\end{table}

\paragraph*{Constant borders}\label{par:borders} The \giraf package allows lattices to have fixed borders. More precisely, each site on the boundary of the lattice has neighbors with respect to the topology induced by $\G$ that are set to constants. Setting the borders can be done using the optional inputs $\textsf{top}=\textsf{c}(t_1,\ldots,t_w)$, $\textsf{left}=\textsf{c}(\ell_1,\ldots,\ell_h)$, $\textsf{bottom}=\textsf{c}(b_1,\ldots,b_w)$, $\textsf{right}=\textsf{c}(r_1,\ldots, r_h)$ and $\textsf{corner}=\textsf{c}(c_1,\ldots, c_4)$ (see Figure \ref{fig:neigh-border}). In the absence of a specified border, these will take the default value \textsc{NULL}.


\section{Recursive computing and partition function}
\label{sec:recursion}

The inherent difficulty of Gibbs random fields arises from the intractable normalizing constant, called the partition function, defined by
\begin{equation}
Z(\psi, \G)  = \sum_{\x\in \X} \exp\left\lbrace - H\left(\x\suchthat\psi,\G\right) \right\rbrace.
\label{eqn:partition}
\end{equation}
The latter is a summation over the numerous possible realizations of the random field $\Xv$, that cannot be computed directly except for small grids and small number of colors $K$.

\giraf contains generalised recursion proposed by \cite{reeves2004} for computing the normalising constant of general factorisable models such as the autologistic or the Potts model (\nameref{ncmrf}). This method, based on an algebraic simplification due to the reduction in dependence arising from the Markov property, applies to lattices with a small number $h$ of rows. Note that the complexity of the algorithm is exponential in the number of rows ($\gO(K^h)$) and linear in the number of columns. Hence the current limitation is $K^h \leq 2^{25}$ for computing the normalising constant.

\begin{figure}[t]
\centering
\begin{minipage}[t]{7cm}
\centering
\input{g4border.tex}\\
(a)
%\caption{Graph $\G_4$}
\end{minipage}
\begin{minipage}[t]{7cm}
\centering
\input{g8border.tex}\\
(b)
%\caption{Graph $\G_8$}
\end{minipage}
\caption{\small First and second order neighborhood graphs $\G$. (a) The four closest neighbors graph $\G_4$. Neighbors of the vertex in black are represented by vertices in gray. (b) The eight closest neighbors graph $\G_8$. Neighbors of the vertex in black are represented by vertices in gray.}
\label{fig:neigh-border}
\end{figure}


\subsection{\textsf{NC.mrf}}
\label{ncmrf}

The function \textsf{NC.mrf} computes the partition function \eqref{eqn:partition} of a $K$-state Potts model defined on a rectangular $h\times w$ lattices ($h\leq w$), with either a first order or a second order dependency structure (see Figure \ref{fig:neigh}) and a small number of rows (up to 25 for 2-states models).

\begin{usage}
\begin{multline*}
\textsf{NC.mrf}(\textsf{h, w, param, ncolors = 2, nei = 4, pot = NULL, top = NULL,} \\
        \textsf{left = NULL, bottom = NULL, right = NULL, corner = NULL})
\end{multline*}
\end{usage}

\begin{inputs} see Section \ref{subsec:inputs}.
\end{inputs}

\begin{example}
<<eval = FALSE>>=
# NC.mrf
# Dimension of the lattice
height <- 8
width <- 10

# Interaction parameter
Beta <- 0.6 # Isotropic configuration
# Beta <- c(0.6, 0.6) # Anisotropic configuration for a first 
                      # order dependency structure (nei = 4).
# Beta <- c(0.6, 0.6, 0.6, 0.6) # Anisotropic configuration for a second
                                # order dependency structure (nei = 8).

# Number of colors. Automatically set to 2 if not specified.
K <- 2 
# Number of neighbors. Automatically set to 4 if not specified.
G <- 4

# Optional potential on sites. Automatically set to NULL if not specified
potential <- runif(K,-1,1)
# Optional borders. Automatically set to NULL if not specified
Top <- Bottom <- sample(0:(K-1), width, replace = TRUE)
Left <- Right <- sample(0:(K-1), height, replace = TRUE)
Corner <- sample(0:(K-1), 4, replace = TRUE)

# Partition function for the default setting
NC.mrf(h = height, w = width, param = Beta)
# When specifying the number of colors and neighbors
NC.mrf(h = height, w = width, ncolors = K, nei = G, param = Beta)
# When specifying an optional potential on sites
NC.mrf(h = height, w = width, ncolors = K, nei = G, param = Beta, 
       pot = potential)
# When specifying possible borders. The users will omit to mention all
# the non-existing borders
NC.mrf(h = height, w = width, ncolors = K, nei = G, param = Beta, 
       top = Top, left = Left, bottom = Bottom, right = Right, 
       corner = Corner)
@
\end{example}


\section{Sampling from a Gibbs random field}
\label{sec:sampler}

Sampling from a Gibbs distribution can be a daunting task due to the correlation structure on a high dimensional space, and standard Monte Carlo methods are impracticable except for very specific cases. For rectangular lattices with a small number of row, \cite{friel2007} have extended the recursive algorithm of \cite{reeves2004} to allow exact sampling from the likelihood (\nameref{exactmrf}). But for larger lattices, Markov chain Monte Carlo (MCMC) methods have played a dominant role in dealing with such intractable likelihood, the idea being to generate a Markov chain whose stationary distribution is the distribution of interest. The \giraf package contains two widely used MCMC procedure to give approximate samples: the Gibbs sampler and the Swendsen-Wang algorithm (\nameref{samplermrf}).

\subsection{\textsf{exact.mrf}}
\label{exactmrf}

The function \textsf{exact.mrf} gives exact sample from the likelihood \eqref{eqn:gibbs} of a $K$-state Potts model defined on a rectangular $h\times w$ lattice ($h\leq w$), with either a first order or a second order dependency structure (see Figure \ref{fig:neigh}) and a small number of rows (up to 19 for 2-states models). More generaly, due to the complexity of the algorithm, the current limitation for exact sampling is $K^h \leq 2^{19}$.


\begin{usage}
\begin{align*}
\textsf{exact.mrf}(\textsf{h, w, param, ncolors = 2,} & \textsf{ nei = 4, pot = NULL,} \\
    \textsf{ top = NULL,} & \textsf{ left = NULL, bottom = NULL,} \\
        &\textsf{ right = NULL, corner = NULL, view = FALSE})
\end{align*}
\end{usage}

\begin{inputs} see also Section \ref{subsec:inputs}.
\begin{center}
\begin{tabular}{p{0.20\textwidth}p{0.74\textwidth}}
\textsf{\bf view} & Logical value indicating whether the draw should be printed. Do not display the optional borders.
\end{tabular}
\end{center}
\end{inputs}

\begin{example}
<<eval = FALSE>>=
# exact.mrf
# Dimension of the lattice
height <- 8
width <- 10

# Interaction parameter
Beta <- 0.6 # Isotropic configuration
# Beta <- c(0.6, 0.6) # Anisotropic configuration for a first 
                      # order dependency structure (nei = 4).
# Beta <- c(0.6, 0.6, 0.6, 0.6) # Anisotropic configuration for a second
                                # order dependency structure (nei = 8).

# Number of colors. Automatically set to 2 if not specified.
K <- 2 
# Number of neighbors. Automatically set to 4 if not specified.
G <- 4

# Optional potential on sites. Automatically set to NULL if not specified
potential <- runif(K,-1,1)
# Optional borders. Automatically set to NULL if not specified
Top <- Bottom <- sample(0:(K-1), width, replace = TRUE)
Left <- Right <- sample(0:(K-1), height, replace = TRUE)
Corner <- sample(0:(K-1), 4, replace = TRUE)

# Exact sampling for the default setting
exact.mrf(h = height, w = width, param = Beta, view = TRUE)
# When specifying the number of colors and neighbors
exact.mrf(h = height, w = width, ncolors = K, nei = G, param = Beta, 
          view = TRUE)
# When specifying an optional potential on sites
exact.mrf(h = height, w = width, ncolors = K, nei = G, param = Beta, 
       pot = potential, view = TRUE)
# When specifying possible borders. The users will omit to mention all
# the non-existing borders
exact.mrf(h = height, w = width, ncolors = K, nei = G, param = Beta, 
          top = Top, left = Left, bottom = Bottom, 
          right = Right, corner = Corner, view = TRUE)
@
\end{example}



\subsection{\textsf{sampler.mrf}}
\label{samplermrf}

The function \textsf{sampler.mrf} gives approximate sample from the likelihood \eqref{eqn:gibbs} of a $K$-state Potts model defined on a rectangular $h\times w$ lattice ($h\leq w$), with either a first order or a second order dependency structure (see Figure \ref{fig:neigh}), using one of the following MCMC samplers.

The Gibbs sampler is a highly popular MCMC algorithm in Bayesian analysis starting with the influential development of \citep{geman1984}. It can be seen as a component-wise Metropolis-Hastings algorithm \citep{metropolis1953, hastings1970} where variables are updated one at a time and for which proposal distributions are the full conditionals themselves.

The Swendsen-Wang algorithm \citep{sw87} originally designed to speed up simulation of Potts model close to the phase transition\footnote{Symmetry breaking for large values of parameter $\beta$ due to a discontinuity of the partition function when the number of sites $n$ tends to infinity. When the parameter is above the critical value $\log\left\lbrace 1+\sqrt{K}\right\rbrace$ of phase transition, one moves gradually to a multi-modal distribution, that is, values $x_i$ are almost all equal} and bypass slow mixing issues of the Gibbs sampler.
This algorithm makes a use of auxiliary variables in order to incorporate simultaneous updates of large homogeneous regions. Note that we solely use Swendsen-Wang updates and we do not run the chain as a \textit{complete coupling} chain as proposed by \cite{huber1999} to get exact sample from the likelihood.




\begin{usage}
\begin{align*}
\textsf{sampler.mrf(iter, sampler = "Gi} & \textsf{bbs", h, w, ncolors = 2, nei = 4, } \\
                     \textsf{param, } & \textsf{ pot = NULL, top = NULL, left = NULL, } \\
                      \textsf{bot}   & \textsf{tom = NULL, right = NULL, corner = NULL,} \\
                        & \textsf{initialise = TRUE, random = TRUE, view = FALSE})
\end{align*}
\end{usage}

\begin{inputs} see also Section \ref{subsec:inputs}.
\begin{center}
\begin{tabular}{p{0.20\textwidth}p{0.74\textwidth}}
\textsf{\bf iter} & Number of iterations of the algorithm.\\
\textsf{\bf sampler} & The method to be used. The latter must be one of "Gibbs" or "SW" corresponding respectively to the Gibbs sampler and the Swendsen-Wang algorithm.\\
\textsf{\bf initialise} & Logical value indicating whether initial guess should be randomly drawn.\\
\textsf{\bf random} & Logical value indicating whether the sites should be updated sequentially or randomdly. Used only with the "Gibbs" option.\\
\textsf{\bf view} & Logical value indicating whether the draw should be printed. Do not display the optional borders.
\end{tabular}
\end{center}
\end{inputs}

\begin{example}
<<eval = FALSE>>=
# sampler.mrf
# Algorithm settings
n <- 200
method <- "Gibbs"

# Dimension of the lattice
height <- width <- 100

# Interaction parameter
Beta <- 0.6 # Isotropic configuration
# Beta <- c(0.6, 0.6) # Anisotropic configuration for a first 
                      # order dependency structure (nei = 4).
# Beta <- c(0.6, 0.6, 0.6, 0.6) # Anisotropic configuration for a second
                                # order dependency structure (nei = 8).

# Number of colors. Automatically set to 2 if not specified.
K <- 2 
# Number of neighbors. Automatically set to 4 if not specified.
G <- 4

# Optional potential on sites. Automatically set to NULL if not specified
potential <- runif(K,-1,1)
# Optional borders. Automatically set to NULL if not specified
Top <- Bottom <- sample(0:(K-1), width, replace = TRUE)
Left <- Right <- sample(0:(K-1), height, replace = TRUE)
Corner <- sample(0:(K-1), 4, replace = TRUE)

# Sampling method for the default setting
sampler.mrf(iter = n, sampler = method, h = height, w = width, 
            param = Beta, view = TRUE)
# Sampling using an existing configuration as starting point
sampler.mrf(iter = n, sampler = method, h = height, w = width, 
            ncolors = K, nei = G, param = Beta, 
            initialise = FALSE, view = TRUE)
# Specifying optional arguments. The users may omit to mention all
# the non-existing borders
sampler.mrf(iter = n, sampler = method, h = height, w = width, 
            ncolors = K, nei = G, param = Beta,
            pot = potential, top = Top, left = Left, bottom = Bottom, 
            right = Right, corner = Corner, view = TRUE)
# Gibbs sampler with sequential updates of the sites. 
sampler.mrf(iter = n, sampler = "Gibbs", h = height, w = width, 
            ncolors = K, nei = G, param = Beta,
            random = FALSE, view = TRUE)
@
\end{example}



\phantomsection
\addcontentsline{toc}{section}{References}
\bibliographystyle{abbrvnat}
\bibliography{bibliography}



\end{document}