\chapter{Introduction}
The last two decades have witnessed remarkable advances in the fields of
biotechnology and computation.  As a result, biologists who study or manage
natural populations of plants and animals today  have access to numerous
new tools.  Of great value are techniques for assaying
genetic variability within individuals, and within populations.  With the
advent of polymerase chain reaction (PCR) and the discovery of new, 
polymorphic, genetic markers, a number of questions about populations and the
interactions of the individuals within them may be addressed using genetic
data.  The advances in computation have been equally significant.  Today's
computers are fast enough to allow numerically intensive statistical analyses
to be run on desktop machines.  This has led to significant development in the
field of computational statistics, particularly in the use of Monte Carlo and
Markov chain Monte Carlo (MCMC) techniques for computing likelihoods and
posterior probabilities.   This dissertation focuses on developing
likelihood-based statistical procedures and the computational methods
necessary to  apply them  to the analysis of genetic data in the context of
questions of interest to biologists, conservation biologists and wildlife
managers.


Each of Chapters~\ref{chap:ForBack} through~\ref{chap:Hybrids} provides an almost
self-contained account of the application of statistical theory and computational
techniques to a practical problem in population genetics.  Chapter~\ref{chap:ForBack}
describes a hidden Markov model and an importance-sampling method for evaluating the
likelihood for the parameter
$N_e$---the genetically effective size---of a population from temporally-spaced
genetic data.  Chapter~\ref{chap:Urns} considers the direct estimation of the ratio
of the effective number of breeders to the census number of breeders; a parameter we
call
$\lambda$.  This chapter describes problems with modeling genetic
transmission in populations via the well-known Wright-Fisher model, and
advances, instead, a model of genetic inheritance derived from a P\'{o}lya urn
scheme.  Using a hidden Markov model with transitions based on this urn scheme,
the final section of the chapter describes an MCMC method for computing the
likelihood or posterior distribution for $\lambda$ or $N_e$. 
Chapter~\ref{chap:overlap} uses the urn scheme of the previous chapter to develop a
likelihood model for
$\lambda$ in the case of salmon populations which, because of their interesting life
histories, violate the discrete generation assumption of the standard Wright-Fisher
model.   This likelihood is used for Bayesian inference of $\lambda$ in salmon
populations.  Once again, this relies on MCMC. 


 Chapters~\ref{chap:MixAdmix} and~\ref{chap:Hybrids}  deal with a different
problem in population  genetics---that of identifying the origin of 
individuals in  mixed and admixed (\ie including some hybrid individuals)
populations.  Chapter~\ref{chap:MixAdmix} describes an extension to the work of
\citeN{Pritchardetal2000} allowing individuals in an admixed population to be
of either pure or admixed genetic origin and uses an example with data from
wildcats ({\em Felis sylvestris}) in Scotland.  Chapter~\ref{chap:Hybrids} develops
a model for admixture that is more appropriate when hybrids have only been formed
over a small number of generations.  This method uses genetic data to compute the
posterior probability that an individual in a sample is, for example, a pure
individual, an $\mathrm{F}_1$ or $\mathrm{F}_2$ hybrid, or a
backcross.  The method is demonstrated on data from coastal cutthroat trout
({\em Oncorhynchus clarki}), steelhead trout ({\em O.~mykiss}), and their
hybrids.   The final chapter, 7, includes brief concluding remarks and
describes  further work to be done in the areas treated by this
dissertation.


This leaves the present chapter to briefly summarize how and why the main
chapters of this dissertation are related and connected.  I see
three levels of connection.   First,  each practical
application described confronts a problem in the
study and management of natural (non-human) populations using genetic data that
are increasingly available.  In my case, each of these applications has been
motivated by a particular problem which occurs in the study of salmon populations. 
Second, each chapter makes heavy use of Monte Carlo computational techniques.  While
Monte Carlo likelihood  and MCMC have been used since the early nineties in
genetic problems involving computations on pedigrees, and soon thereafter were
applied to inference on evolutionary time scales, their application to
population genetics problems has been more recent, and the chapters herein
represent some of the first applications of MCMC to genetic inference at the
population level.  The third point is more of a statistical curio than a
unifying theme, but I find it interesting nonetheless---despite the
fact that the problems addressed in Chapters~\ref{chap:MixAdmix}
and~\ref{chap:Hybrids} are quite different from those in Chapters~\ref{chap:ForBack}
through~\ref{chap:overlap}, we encounter in all of them useful applications of hidden
Markov models, mixture models, and urn models.  In the remaining sections of this
chapter, I will elaborate on the first two points, and for those unfamiliar with
them, I provide a brief introduction to Monte Carlo techniques.

\section{Genetic Data}
Advances in biotechnology have revolutionized conservation biology and resource
management.
\citeN{Aviseetal1995} review the genetic markers currently
available to researchers, discuss the types of analyses those
markers allow, and review applications in conservation genetics.  
The advent of  Mendelian-inherited microsatellite markers
\cite{Tautz1989,Wri&Ben1994} has made informative genetic data
increasingly available and inexpensive for such applications. 
Among other examples, DNA markers amplified from fin clips
have been used in monitoring Pacific salmon
\cite{Olsenetal1996}, while hair samples have
been used in studying bear
\cite{Taberletetal1997}, 
%%pandas \cite{Zhangetal1994}, macaques\cite{Ino&Tak1993}, 
and chimpanzee \cite{Morinetal1993} populations.  
PCR-based technologies are especially appropriate for populations of
conservation interest as sampling is non-destructive and/or
non-invasive. It is
thus possible to obtain data at multiple time points, and,
since the markers are typically polymorphic, the data are
informative in  characterizing the population at each time
point, and hence also in detecting and quantifying the gene
frequency changes  caused by small effective population size or
genetic exchange with other populations.


With microsatellite markers and PCR, data may be extracted from 
archived tissues, giving the opportunity to
obtain  data from time points in a population's
past.  For example, museum-preserved skins from known
populations of the pocket gopher 
provide genetic data on the
populations at two time points (1950's, 1970's) which may be
compared to current samples (Ellie Steinberg, UW Dept.\ of
Zoology, pers.\ comm.).  For some fish populations, the 
situation is even better.   Some such populations have been the
subject of long-term ecological research efforts with
population size estimates available on a yearly basis, and age
composition inferred from fish scales.  Genetic marker data may
be obtained from these archived scales.  Recently,  
\citeN{Mil&Kap1997} isolated DNA from northen pike scales
collected from Lake Escanaba, WI. Using data from three
years, 1961, `77, and `93, they estimated $N_e$ from the
temporal changes in allele frequencies over the two time
intervals.  In a similar, ongoing study, microsatellites from
archived juvenile  Keogh River (Vancouver Island) and Snow Creek
(Washington State) steelhead scales provided ample material for amplifying
microsatellite markers by PCR \cite{Ardren1999}.  The use of archived fish tissues
for population genetic studies is further described by \citeN{Nielsenetal1999}.
 


\section{Salmon Populations and Genetics}
Pacific salmon ({\em Oncorhynchus} spp.)\ and their relatives, 
Atlantic salmon and the true trouts ({\em Salmo} spp.)\ and chars ({\em Salvelinus}
spp.)\ in the family Salmonidae are an evolutionarily fascinating and commercially
valuable group of fish.  Both of these factors have contributed to the generation of
an enormous amount of genetic data on these species.  The five species of
Pacific salmon native to the West Coast of North America, sockeye ({\em O.~nerka}),
chinook ({\em O.~tshawytscha}), coho ({\em O.~kisutch}), chum ({\em O.~keta}), and
pink ({\em O.~gorbuscha}) exhibit an anadromous and semelparous life history; they
hatch from eggs in fresh water, migrate to the ocean and mature to adulthood,
then return to their natal fresh waters to spawn, and then die soon thereafter.  A
remarkable feature of these migrations is the salmon's homing ability.  Returning
adults typically spawn in the stream in which they
were born.  While this homing behavior is not always perfect, and straying (spawning
in a non-natal stream) is known to occur, their homing does lead to a situation in
which salmon populations represent relatively isolated, reproducing populations with
low gene flow between populations.  A comprehensive review of life histories in
Pacific salmon may be found in the volume edited by \citeN{Gro&Mar1991}.

  Atlantic salmon
({\em Salmo salar}), steelhead trout ({\em O.~mykiss}), and coastal
cutthroat trout ({\em O.~clarki}) display anadromous forms and homing ability;
however, they do not necessarily die after spawning.  {\em O.~mykiss} and {\em O.
clarki} also exhibit non-anadromous or  resident forms.   Many of the chars ({\em
Salvelinus} spp.) also exhibit resident and anadromous forms.  There has been
considerably less research on chars than on Pacific salmon and the true
trouts, but there is still considerable genetic information about them (see  
\citeN{Lea&All1997} and \citeN{Tayloretal2001} and references therein).

Commercial fishing for Pacific salmon is a mammoth industry on the West Coast. Sport
fishing is also economically important to the region.  These influences have had a
siginficant impact on salmon population abundnace.  Additionally, the nature of
their life history makes salmon vulnerable to a number of different
anthropogenic disturbances that affect the environments in which they live and
through which they migrate.  For example, hydropower  operations,
forestry, agriculture, road-building, and wetlands destruction all impact salmon
at some some point in their life history.  In response to dwindling population sizes
and the wholesale destruction of spawning habitat for hydroelectric power
generation, a great number of salmon hatcheries have been built to try to maintain
or supplement salmon production on the West Coast. In recent years the
aquaculture practice of net-pen rearing of Atlantic salmon on
both coasts has increased dramatically.  Hatchery and aquaculture practices have
brought with them their own suite of impacts on wild salmon populations. In the
lower 48 states of the United States, in particular, but also in Canada and Alaska,
all of these impacts have resulted in extinction or crisis for a significant number
of salmon populations \cite{Nehlsenetal1991}.

The use of genetic markers in the management of salmonid stocks has a long and
extensive history
\cite{Rym&Utt1987}.  More recently, with the 
application of the Endangered Species Act to populations of Pacific salmon, genetic
characteristics of populations along with other traits are used in delineating
the ``Evolutionarily Significant Units" to which the Act is applied
\cite{Waples1995}.  Consequently, there is a vast literature on salmon genetics and
the use of genetic markers in the conservation of salmon.  My purpose in the
remainder of this section is not to attempt a comprehensive review of that
literature (such reviews and case studies may be found in \citeN{All&Wap1996} and
the articles co-published with \citeN{Utter1999}) but rather to first give the
reader some sense of the varieties of genetic data on salmon and then to
highlight the uses of those data that are particularly relevant to this
dissertation.  

The first, and still the most widely-used, genetic markers for
salmon are electrophoretically detectable enzymes described and
reviewed by \citeN{Utteretal1987}.  Today there are well over 60 commonly-used
allozyme loci routinely used by the National Marine Fisheries Service in screening
salmon populations for genetic variation and for other purposes.  However, many 
newer, molecular markers have been discovered and used in salmon as well. 
\citeN{Gyl&Wil1987} describe mitochondrial DNA (mtDNA) markers in the era before
PCR.  Also in the pre-PCR era, various probes were developed for hybridization to
restriction-enzyme-digested total genomic DNA of salmon 
\cite{Devlinetal1991}.  With the ability to amplify DNA by PCR, mtDNA markers were
more easily prepared and other, new types markers became available. First were
markers that were associated with known, functional genes like those coding for
ribosomal DNA
\cite{Pendasetal1995} or growth hormone
\cite{Gro&Nil1995}.  Soon thereafter, however, probes detected minisatellite
loci in all salmon species \cite{Prodohletal1995}, and the method of Random
Amplification of Polymorphic DNA (RAPD) was applied to salmon \cite{Eloetal1997}. 
Today, microsatellite markers \cite{Olsenetal1996} are widely available for salmon. 
Research into new markers for salmon continues, with the discovery of new short
interspersed repeat (SINE) segments being a recent example \cite{Perezetal1999}.

Genetic markers have been used in at least several hundreds of studies of salmon
populations.   Some of these concern themselves primarily with reporting genetic
variation across and within populations; however, many of them are directed toward
answering specific questions or estimating particular quantities associated with the
populations.  Of particular relevance to this dissertation are studies involving
genetic stock identification, the estimation of effective population size, and the
detection of hybrids.  In  later chapters, I present novel computational
methods for each these tasks.  

Genetic stock identification refers to the use of multilocus genotypes (without
knowledge of phase) to assign individuals sampled from a mixture of fish
to one of the populations contributing to the mixture, and also to estimate the
proportion of fish in the mixture from each of the source populations.   The
empirical studies employing these techniques are too numerous to list; however,
\citeN{Woodetal1989} is a representative example using allozymes and other
biological trait data, \citeN{Beacham1996} describes genetic stock identification
using minisatellites  and
\citeN{Beachametal2000} and
\citeN{Olsenetal2000} both provide representative cases using microsatellite
markers.  


Interestingly, the  availability  in the early 1980's of allozyme
data from salmon populations and problems in salmon population management were the
primary factors motivating the original development of statistical methodologies for
these sorts of genetic mixtures. \citeN{Milneretal1981} present the basic
framework and an EM-algorithm for maximum likelihood methods in genetic stock
identification conditional on allele frequencies being known from the populations
contributing to the mixture.  \citeN{Millar1987} analyzes conditions of
identifiability and shows that another method known as the ``classification" method
is a special case of the the likelihood approach.  \citeN{Smouseetal1990}
present an EM-algorithm for inference in the case when the allele frequencies from
the contributing stocks is not assumed known without error, and \citeN{Pel&Mas2001}
give a Bayesian version using MCMC to simulate from the posterior distributions
of interest.

The estimation of effective size of salmon populations is another problem that has
spurred the development and refinement of statistical methods.  The work of
\citeN{Waples1989}, clarifying and generalizing previously developed $F$-statistic
methods for estimating effective population size from temporal changes in allele
frequencies, was motivated in large part by problems in salmon biology.  In a
series of papers, he explores genetic change over time
\cite{Wap&Tee1990,Waples1990(2)} in salmon populations and tailors an
$F$-statistic method for estimating effective size to the Pacific salmon
life history \cite{Waples1990}. This method was later made more practical by an
algorithm presented in
\citeN{Tajima1992}.  These methods have been used to estimate the effective sizes of
endangered salmon populations
\cite{Waplesetal1993}.  Estimation of brown trout ({\em Salmo trutta})  effective
sizes \cite{Jor&Rym1996} motivated the development of another method for applying
$F$-statistic  estimators to temporal allele frequency data for the estimation of
effective size in populations with overlapping generations \cite{Jor&Rym1996}.  

A number of empirical studies have been conducted, estimating the effective size or
the effective number of spawners in salmon populations.  \citeN{Ardren1999} uses
microsatellites amplified from fish scales to try to estimate the effective size of
two West Coast steelhead populations.  \citeN{Kincaid1995} analyzes mating patterns
and breeding history to try to estimate inbreeding effective size in hatchery
populations of salmonids.  \citeN{Hedricketal1995} discuss and estimate the
effective size of winter run chinook salmon.  \citeN{Tes&Ber1997}, using
microsatellite and mtDNA markers, assess the effective population size of landlocked
Atlantic salmon populations and the effect that supportive breeding programs have
on the effective size of these populations.  

The literature on the use of genetic methods for detecting hybrids between salmonid
species is enormous and will not be exhaustively reviewed here.  The studies using
genetic markers include those that investigate naturally-occurring  
hybridization between sympatric species (\cf
\citeNP{Cam&Utt1985}; \citeNP{Tayloretal2001}; \citeNP{Eloetal1995}) and numerous
ones involving interbreeding between hatchery or farm-reared salmon and
native populations of conspecifics (for example, \citeNP{Cliffordetal1998}) or other
species
\cite{Jan&Oes1997}.

\section{Monte Carlo in Genetics}
Many inference problems in statistical genetics involve complex stochastic models
that include a great number of variables.  Typically only a fraction of these
variables can be directly observed.  These variables summarize the observed genetic
data.  The remainder of the variables in the stochastic model are not directly
observable and are referred to as {\em latent} variables.  The likelihood function
for such inference problems  can be expressed as the sum over the latent
variables of the joint probability  of the observed data and  the latent variables,
conditional on the genetic parameters of interest.  Often, however, the space of
latent variables is huge and that sum is not directly computable. 

Monte Carlo methods are stochastic integration techniques that are useful
for approximating such intractable sums.  I provide a brief introduction to Monte
Carlo methods in Section~\ref{sec:MCIntro}, after establishing conventions for
mathematical notation.  Here I briefly review the recent history of Monte Carlo
techniques, and, in particular, Markov chain Monte Carlo (MCMC) techniques in
statistical genetics.  

Four years ago, as I was starting the research reported in this dissertation, 
MCMC methods had been used extensively  
in the analysis of data on extended or complex pedigrees.
In segregation and linkage analysis of trait and genetic marker data 
observed on members of a known pedigree, MCMC is used to sample over the
space of latent variables which may be defined as genotypes
\cite{Guo&Tho1992} or meiosis indicators
\cite{Thompson1994} or both \cite{Lan&Mat1989}.  Much work had been done
in this field to improve the mixing of chains in the space of latent variables.
One example in the area of inference of ancestral types on a large complex
pedigree structure is the simulated tempering method of \cite{Gey&Tho1995}.
To improve mixing and ensure irreducibility, methods were developed in which
multiple components of the latent variables are updated  simultaneously.  In pedigree
analysis, such methods include use of a block-updating Gibbs sampler
\cite{Janssetal1995}, a whole-meiosis Gibbs sampler \cite{Tho&Hea1999}, and a
whole-locus Gibbs sampler
\cite{Kong1991}.  
%% These methods enable efficient computations on extended pedigrees.  
Brand new MCMC methods, at the time, like  reversible-jump MCMC
\cite{Green1995} were quickly adopted to analyze more
complex  model spaces in genetic analysis.  \citeN{Heath1997} applies reversible jump
methods to detect and locate  multiple quantitative trait loci (QTL) 
from trait and genome-scan data, where the number of QTL is not 
prespecified and thus the dimension of the model varies within a single MCMC run.

%% Thus the dimension of the genetic model, which depends on 
%% the number of QTL, may vary within a
%% single run of the Markov chain 
%% (\citeNP{Satagopanetal1996}, \citeNP{Heath1997}).
 
MCMC methods had also been used in analyses of inference of relationship 
among individuals from genetic data.  
\citeN{Painter1997} develops methods for estimation 
of sibship structure, sampling directly over the space of alternative 
sibship structures, using data on microsatellite markers.
\citeN{Geyeretal1993} use a Metropolis-Hastings MCMC 
method to construct a Monte Carlo likelihood function for relationship 
parameters among a group of individual California condors,
on whom there are multilocus DNA fingerprint data.  At the other extreme
of the evolutionary time scale, 
MCMC methods had also been used in phylogenetic analyses, to estimate 
evolutionary parameters, such as the product of mutation rate and 
effective size
\cite{Kuhneretal1995}, or the rate of increase of populations
\cite{Kuhneretal1997}.
In these analyses, the latent variable is the structure 
and inter-coalescence times of the ancestral coalescent
\cite{Kingman1982}
of a sample of DNA sequences, and is sampled
using a Metropolis-Hastings algorithm. 
\citeN{Newtonetal1997} propose an alternative specification of the 
coalescent structure, leading to an MCMC sampler which can make large 
changes in ancestral topology
in a single MCMC step. In some cases, this specification may
provide a better mixing sampler.
Other Monte Carlo likelihood methods
had also been used for coalescent models in the context of estimation of
growth rates \cite{Gri&Tav1994(2)} and recombination rates
\cite{Gri&Maj1996},  and inference of mutation models
\cite{Nielsen1997}.  


%Again the realized latent variable is the
%ancestral coalescent, in these cases specified by the 
%sequence of ancestral recombination, mutation, or coalescent events. 
%\citeN{Ste&Don2000} recently derived a far superior importance
%sampling method in such problems.

Since that time, the use of MCMC methods in statistical genetics has grown
dramatically, due in large part to the growing acceptance of Bayesian techniques and
the close association between Bayesian computation and MCMC\@.  \citeN{Stephens2001}
reviews recent MCMC and importance sampling methods for likelihood and Bayesian
inference using coalescent models, and the numerous new MCMC approaches for
inference of phylogenies are reviewed by \citeN{Hue&Bol2001}.  MCMC methods have also
enjoyed further refinements and applications in the detection of quantitative
trait loci from data on outbred pedigrees \cite{Hoeschele2001}  and in the inference
of breeding values in animal breeding \cite{Gianola2001}. 



While MCMC methods have been used on pedigrees and on coalescents for some time,
only recently have they been used in conservation genetics problems at the
intervening population time-scale where a pedigree structure is not available and
the data, either observed or latent,  are allele frequencies in specified
populations.  The methods developed in this thesis are applications of MCMC to
problems relevant in conservation genetics.  Other Monte Carlo approaches relevant to
conservation genetics on short time scales have been developed in the last several
years and include the detection and characterization of recent bottlenecks using
genetic data \cite{Beaumont1999}, the estimation of effective population size
\cite{Kitadaetal2000}, the analysis of population structure and admixture
\cite{Pritchardetal2000}, and genetic stock identification in fisheries management
\cite{Pel&Mas2001}.


\section{Notational Conventions}
Statistical genetics is rich with complex data arrangements and, therefore,
multiply subscripted variables.  To avoid confusion, I adopt the following
conventions for mathematical notation, and apply them throughout the chapters
of the dissertation:  scalar random variables are either uppercase, italicized, Roman
letters or lowercase Greek letters; for example $X$ or $\theta$. Greek letters are
also used for quantities considered to be parameters in a frequentist setting.   For
the Roman letters, vector-valued random variables are written in slant-bold, with
their components being written as scalars with a single subscript:
$\btX = (X_1,\ldots,X_n)$.  Random variables which are collections of vectors are
denoted by Roman bold, $\bX$, with the
$i\thh$ vector element within it denoted  by $\btX_i$ and scalar components by
double subscripts separated by a comma, \ie $\bX = (\btX_1,\ldots,\btX_m)$
with $\btX_i = (X_{i,1},\ldots,X_{i,n})$.  Some collections of
vectors require more than one subscript to denote each vector.  These
are still denoted by Roman bold characters with the obvious extensions; for
example  $\bX = (\bX_1,\ldots,\bX_k)$, $\bX_t = (\btX_{t,1},\ldots,\btX_{t,m})$,
and $\btX_{t,j} = (X_{t,j,1},\ldots,X_{t,j,n})$.  The shading and
subscripting conventions for Greek letter random variables or parameters are
similar, with uppercase bold substituted for collections of vectors:  scalar
$\theta$, vector
$\bttheta$, collection of vectors (and beyond) $\bTheta$.
 

$P(\cdot)$  denotes probability mass and density
functions alike.  $P(X)$ is used as a name for the distribution
of the variable $X$ and also to denote the marginal
probability mass or density at a realized value of $X$, \ie 
$P(X)$ also serves as a shorthand for $P(X=x)$.  Only in cases where
confusion would be likely or where the distinction is of central importance shall I
make a notational distinction between the realized value of a random variable and
the random variable itself.  In such cases, as in the following section on Monte
Carlo, I use the lowercase version of the random variable to denote the realized
value.  When dealing with a sequence of realized values of a random variable $X$
generated as a sample to be used for Monte Carlo, I denote the $i\thh$ such member
of the sample by
$X^{(i)}$.  The notation $P(\cdot|\cdot)$ is used for conditional probability mass
and density functions.  Hence $P(X|Y)$ is the conditional probability of $X$ given
$Y$. When multiple variables are involved, commas are used between them;
$P(Y,Z|W,X)$ is the conditional joint probability of $Y$ and $Z$ given $W$ and $X$.
Sometimes the probability depends on a parameter, and that dependence is more
conveniently expressed by subscripting the
$P$.  For example,
$P_\theta(X|Y)$.   The expected value of a random variable $X$, is written
$\Exp(X)$, and the conditional expectation of $X$ given $Y$ is written $\Exp(X|Y)$. 
If the expection of $X$ is taken with reference to a distribution indexed by a
particular parameter, say
$\theta_0$, it may be written as
$\Exp_{\theta_0}$.  The variance of $X$ is denoted $\Var(X)$.

Occasionally, such as in the following section or while denoting proposal
distributions for Metropolis Hastings sampling,  I shall use upper or lowercase
italic letters to refer to probability mass or density functions of
random variables.  


\section{Monte Carlo techniques\label{sec:MCIntro}}
This section provides  elementary background on the Monte Carlo method,
importance sampling, Rao-Blackwellization, and Markov chain Monte Carlo.

\subsection{History and definition of Monte Carlo} The term ``Monte Carlo" was
apparently used by Stanislaw Ulam and John von Neumann as a Los Alamos code
word for the stochastic simulations they applied to building hydrogen
bombs following World War~II.  Shortly after the War, and coinciding with the debut
of the ENIAC computer in 1947, von Neumann and Ulam suggested that the ENIAC would
be useful for applying ``statistical sampling" approximations to solving the
problem of neutron diffusion in fissionable material.  Their methods, involving the
laws of chance, performed well, and soon thereafter were aptly named by Nick
Metropolis after Monte Carlo, the international gaming destination. The moniker stuck
and soon after the War a wide range of difficult problems yielded to the new
techniques \cite{Metropolis1987}.  Despite the widespread use of the methods, and
numerous descriptions of them in articles and monographs, it is virtually impossible
to find a succint definition of ``Monte Carlo method" in the literature.  Perhaps
this is owing to the intuitive nature of the topic which spawns many definitions by 
way of specific examples.  Some authors prefer to use the term ``stochastic
simulation" for almost everything, reserving ``Monte Carlo" only for Monte Carlo
integration and Monte Carlo tests
\cite{Ripley1987}.  Others seem less concerned about blurring the
distinction between simulation studies and Monte Carlo methods.  
Be that as it may, I adopt the following terse definition:
\begin{quotation}
\noindent Monte Carlo is, in essence, the approximation of an
expectation  by the sample mean of
a function of simulated random variables.
\end{quotation}
We will find that this definition is broad enough to cover everything
that has been called Monte Carlo, and yet makes clear the fundamental feature of
Monte Carlo in very familiar terms: Monte Carlo is about invoking laws of large
numbers to approximate expectations.   This applies when the simulated variables
are independent of one another, and may apply when they are correlated with one
another, for example if they are states visited by an ergodic Markov chain.  The
Monte Carlo method is useful precisely because very many quantities of interest may
be expressed as expectations.

While most Monte Carlo simulations are done by computer today, there
were several applications of Monte Carlo methods using coin-flipping,
card-drawing, or needle-tossing (rather than computer-generated
pseudo-random numbers) centuries ago---long before
the name Monte Carlo arose.  

\paragraph*{In more mathematical terms:} Consider a  random variable $X$ (though
depicted as a scalar, all of the following extends to multidimensional random
variables) having probability mass function or probability density function
$f_X(x)$ which is greater than zero on a set of values
$\mathcal{X}$.  Then the expected value of a function $g$ of $X$ is
\begin{equation}
\Exp[g(X)] = \sum_{x\in\mathcal{X}} g(x)f_X(x)
\end{equation}
if $X$ is discrete, and
\begin{equation}
\Exp[g(X)] = \int_{x\in\mathcal{X}} g(x)f_X(x)dx
\end{equation}
if $X$ is continuous.  If we were to take an $n$-sample of  $X$'s,
 $(x^{(1)},\ldots,x^{(n)})$\footnote{In this section, the distinction between the
random variable $X$ and its realized value $x$ is crucial, so a notational
distinction is made between the two---capital $X$ and $X^{(i)}$ refer to random
variables while lowercase $x$ and $x^{(i)}$ refer to realized values of the
random variable $X$.},   and we computed the mean of
$g(x)$ over the sample, then we would have the Monte Carlo {\em estimate} 
\[
\meang =
\frac{1}{n}\sum_{i=1}^n g(x^{(i)})
\]
 of $\Exp[g(X)]$. We could, alternatively, 
speak of the random variable 
\[
\Meang = \frac{1}{n}\sum_{i=1}^n g(X^{(i)}) 
\]
which we call the Monte Carlo {\em estimator} of $\Exp[g(X)]$.


If $\Exp[g(X)]$, exists, then the weak law of large
numbers tells us that for any arbitrarily small $\epsilon$
\[
\lim_{n\rightarrow\infty}P\left(\left|\Meang -
\Exp[g(X)]\right|\ge\epsilon\right) = 0.
\]
This tells us that as $n$ gets large,
there is small probability that $\Meang$ deviates more than a tiny bit from
$\Exp[g(X)]$.  For our purposes, the strong law of large numbers says
much the same thing---the important part being that so long as $n$ is
large enough, $\meang$ arising from a Monte Carlo experiment shall be
as close to $\Exp[g(X)]$ as desired.  This
extends to samples from Markov chains via the weak law of large numbers
for the number of passages through a recurrent state in an ergodic Markov
chain (see \citeNP{Feller1957}).

\pagestyle{myheadings}

It should also be clear  that
$\Meang$ is unbiased for  $\Exp[g(X)]$:
\[
\Exp[\Meang] = \Exp\left(\frac{1}{n}\sum_{i=1}^n g(X^{(i)})\right) =
\frac{1}{n} \sum_{i=1}^n \Exp[g(X^{(i)})] = \Exp[g(X)].
\]

These properties of random samples become useful when one realizes that very many
quantities of interest may be cast as expectations.  Most importantly for
applications in this dissertation, it is possible to express all probabilities and
summations as expectations:
\begin{description}
\item[Probabilities:]  Let $Y$ be a random variable. The probability
that $Y$ takes on some value in a set $\mathcal{A}$ can be expressed as an
expectation using the indicator function:
\begin{equation}
P(Y\in \mathcal{A}) = \Exp[\ind\{Y\in\mathcal{A}\}]
\label{eq:ProbMC}
\end{equation}
where $\ind\{Y\in\mathcal{A}\}$ is the indicator function that takes the value 1
when
$Y\in \mathcal{A}$ and 0 when $Y\not\in \mathcal{A}$.
\item[Summations:] Any sum, even the sum of a collection of 
deterministic variables may be represented as an expectation.  
For example the sum of a function $q(b)$ over the values of $b$
in a finite set
$\mathcal{B}$, may be written as an expectation of a random variable, say $W$,
which takes values in 
$\mathcal{B}$.  In the simplest form, $W$ could take any value in $\mathcal{B}$
with  equal probability 
$p$, and the sum could be cast as the
expectation 
\[
\sum_{b\in B} q(b) = \frac{1}{p}\sum_{b\in B} q(b) p =
\frac{1}{p}\Exp[q(W)].
\]
\end{description}
The immediate consequence of this is that all probabilities
and summations can be approximated by the Monte Carlo method. And further, there is
 no restriction that says
$W$ above must have a uniform distribution.  This was just for easy
illustration.  We will explore this point more while
considering importance sampling.

The variety of quantities that may be estimated by Monte Carlo is great.  In
addition to probabilities and summations, it is also possible to estimate
integrals, probability distributions, and variances, \etc  In each of these cases,
the fundamental feature is the same---the quantity of interest may be expressed as an
expectation which is then approximated by Monte Carlo.  For example, to
approximate a probability distribution, (\ref{eq:ProbMC}) may be applied $K$ times to
$K$ different sets $\mathcal{A}_k,~k=1,\ldots,K$, which, taken together, give a
histogram representation of the distribution.  Or, to estimate the variance of a
random variable $Y$, that variance may be expressed as an expectation (namely the
expected value of $(Y-\Exp Y)^2$) and estimated via Monte Carlo accordingly.  

Many problems in statistical genetics provide examples where a
quantity of interest is a summation.  In such cases the probability
$P(Y)$ of an observed event $Y$ must be computed as the sum over very many latent
variables $X$ of the joint probability $P(Y,X)$.  In such
a case, $Y$ is typically fixed, \ie we have observed $Y=y$, and we are
interested in $P(Y=y)$, but we can't observe the values of the latent
variables which may take values in the space $\mathcal{X}$.  Though it
follows from the laws of probability that
\[
P(Y=y) = \sum_{x\in \mathcal{X}} P(Y=y,X=x),
\]
quite often $\mathcal{X}$ is such a large space (contains so many
elements) that it is impossible to compute the sum.  Application of the
law of conditional probability, however, gives
\begin{equation}
P(Y=y) = \sum_{x\in \mathcal{X}} P(Y=y,X=x) = \sum_{x\in \mathcal{X}}
P(Y=y|X=x)P(X=x).
\label{eq:StandardLatent}
\end{equation}
The term following the last equals sign is the sum over all $x$ of a
function of $x$ (namely, $P(Y=y|X=x)$), weighted by the marginal
probabilities $P(X=x)$.  Clearly this is an expectation, and therefore
may be approximated by Monte Carlo, giving us
\begin{equation}
P(Y=y) \approx \frac{1}{n} \sum_{i=1}^nP(Y=y|X=x^{(i)})
\label{eq:naiveMC}
\end{equation}
where $x^{(i)}$ is the $i\thh$ realization from the marginal distribution of
$X$.  

Unfortunately, (\ref{eq:naiveMC}) would probably provide a very poor Monte Carlo
estimate. Though it is typically easy to formulate a quantity as an expectation and
to propose  a ``naive" Monte Carlo estimator, it is quite another thing to 
have the Monte Carlo estimator  actually provide good estimates in a
reasonable amount of
computer time.  For most problems, a number of Monte Carlo estimators may
be proposed; however some Monte Carlo estimators are clearly better than
others.  Typically, a ``better" Monte Carlo estimator has smaller
variance (for the same amount of computational effort) than its
competitors.  The variance of a Monte Carlo estimator is easily defined.
Going back to our original notation, we have the random variable
$\Meang$, a Monte Carlo estimator of $\Exp(g(X))$.  Like all random
variables, we may compute its variance (if it exists) by the standard
formulas:
\begin{equation}
\Var(\Meang) = \Var\left(\frac{1}{n}\sum_{i=1}^n
g(X^{(i)})\right) = 
\frac{\Var(g(X))}{n} =\frac{1}{n}\sum_{x\in\mathcal{X}}
[g(x) -
\Exp(g(X))]^2f_X(x)
\label{eq:discvar}
\end{equation}
if $X$ is discrete, and
\begin{equation}
\Var(\Meang) = \Var\left(\frac{1}{n}\sum_{i=1}^n
g(X^{(i)})\right) = 
\frac{\Var(g(X))}{n}=\frac{1}{n}\int_{x\in\mathcal{X}} [g(x)-
\Exp(g(X))]^2f_X(x)dx
\label{eq:contvar}
\end{equation}
if $X$ is continuous.  For the rest of this section, we will do everything in
terms of integrals over continuous variables, but it applies equally
well to sums over discrete random variables. 

There are many
ways to reduce the variance of Monte Carlo estimators. Of these
``variance-reduction" techniques, the one
called ``importance sampling" is particularly useful.  I include a short section on
it here, as it will be used in Chapter~\ref{chap:ForBack}. 


\subsection{Importance sampling \label{sec:ImpSampIntro}}
 Importance sampling  \cite{Ham&Han1964} is the art of choosing a good
distribution from which to simulate one's random variables. It involves
multiplying the integrand by 1 (usually dressed up in a ``tricky
fashion") to yield an expectation of a quantity that varies less than the
original integrand over the region of integration.  For example, let
$h(X)$ be a density for the random variable
$X$ which takes values only in $ \mathcal{A}$ so that $\int_{x\in A}h(x)dx = 1$. 
Then
$\frac{h(x)}{h(x)}$ is a ``tricky way" to write 1, and so it follows that 
\begin{equation}
\int_{x\in  \mathcal{A}} g(x)dx = \int_{x\in \mathcal{A} } g(x)\frac{h(x)}{h(x)}dx =
\int_{x\in  \mathcal{A}} \frac{g(x)}{h(x)}h(x) dx =
\Exp_h\left(\frac{g(X)}{h(X)}\right),
\label{eq:ImpSampIdent}
\end{equation}
so long as $h(x)\neq 0$ for any $x\in  \mathcal{A}$ for which $g(x) \neq 0$, and
where
$\Exp_h$ denotes the expectation with respect to the density $h$.  This gives a
Monte Carlo estimator:
\begin{equation}
\widetilde{g_n^h}(X) = \frac{1}{n}\sum_{i=1}^n \frac{g(X^{(i)})}{h(X^{(i)})}~~~~
\mathrm{where}~~~~~X^{(i)} \sim h(X).
\label{eq:impsum}
\end{equation}


Using (\ref{eq:contvar}) and the Cauchy-Schwarz
inequality, it can be shown that $\Var(\gh(X))$ is minimized when $h(x)
\propto |g(x)|$ (see \citeNP{Rubinstein1981}, p.~123).  If we restrict our
attention to what for most of our purposes is the relevant
case,\footnote{This is typically the relevant case because we are interested in
non-negative quantities like probabilities.} that is,
$g(x)\geq 0~\forall x\in  \mathcal{A}$, then it is immediately apparent that the
choice of the density $h(x)$ which minimizes Monte Carlo variance is
proportional to
$g(x)$, \ie if $\alpha h(x) = g(x)$ where $\alpha$ is some constant
of proportionality, then clearly we have $g(x)/h(x)=\alpha~\forall\{x:h(x)>0\}$ so
$\Exp(g(X)/h(X)) =\alpha$ and hence the Monte Carlo variance would be zero by
(\ref{eq:contvar}).

This seems wonderful---to obtain a Monte Carlo estimator with zero
variance we could use (\ref{eq:impsum}), choosing our density $h$ 
proportional to the function $g$.  The absurdity of this wishful
thinking is that the ability to simulate {\em independent} random
variables from $h(x)$, or the ability to compute the density $h(x)$, itself,
implies that the normalizing constant of the distribution is computable, which in
turn would imply that the original integral involving $g(x)$ is computable and
there would hence be no reason to do Monte Carlo at all! Ultimately, however, it
makes clear that a good importance sampling function (as $h$ is called) will be
one that is as close as possible to being proportional to $g(x)$

In summary, a good importance sampling function $h(x)$ has the following
properties:
\begin{enumerate}
\item $h(x) > 0$ whenever $g(x) \neq 0$,~~ (this is required for
(\ref{eq:ImpSampIdent}) to hold)
\item $h(x)$ should be close to being proportional to $|g(x)|$
\item it should be easy to simulate values from $h(x)$
\item it should be easy to compute the density $h(x)$ for any value $x$ that you
might realize.
\end{enumerate}
Fulfilling this wish-list in high dimensional space (where Monte Carlo techniques
are most useful) is often a tall task---it is the main difficulty addressed in
Chapter~\ref{chap:ForBack}.


Note also that $g(x)$ is any arbitrary function, so it certainly includes the
integrand of a standard expectation.  For example, with $X\sim f_X$ we might be
interested in $\Exp(r(X))$ for some function $r$ so we could use 
\[
\Exp(r(X)) = \int r(x)f_X(x)dx = \int\frac{r(x)f_X(x)}{h(x)}h(x) = \Exp_h\left(
\frac{r(x)f_X(x)}{h(x)}\right)
\]
and approximate that by Monte Carlo, simulating values $x^{(1)},\ldots,x^{(n)}$ from
a distribution
$h(x)$ that is close to proportional to
$r(x)f_X(x)$.

Going back to the sum over latent variables problem often
encountered in statistical genetics, importance sampling gives us a
way to improve upon (\ref{eq:naiveMC}). 
From (\ref{eq:StandardLatent}) it is
clear that the optimal importance sampling function would be the conditional
distribution of
$X$ given $Y$, \ie
\[
P(Y=y) = \sum_{x\in \mathcal{X}} P(Y=y,X=x) = \sum_{x\in \mathcal{X}} 
\frac{P(Y=y,X=x)}{P(X|Y=y)}P(X|Y=y).
\]
Note that the right side is a conditional expectation of a function of $X$.
As before $P(X|Y)$ is not computable.  So one must turn to finding some other
distribution, say $P^*(X)$, that is close to $P(X|Y)$ but which is more easily
sampled from and computed.

It is worth noting that, for making likelihood inference of a quantity, it is
sufficient to be able to determine the likelihood function up to an unkown
normalizing constant, \ie it suffices to calculate $\psi P(Y=y)$ where $\psi$ is
unknown.  If the value of the normalizing constant $\psi$ is not required, then the
unnormalized probability
$\psi P(Y=y)$ may be obtained by importance sampling using an importance sampling
function $h(x)$ which is also known only up to a normalizing constant. 
Estimating likelihood ratios in this way has been described by
\citeN{Tho&Guo1991}.    

\paragraph{A common pitfall of importance sampling:}

As a final word on importance sampling, it should be pointed out that {\em the
tails of the distributions matter}\/!  While $h(x)$ might be roughly the same
shape as
$g(x)$, serious difficulties arise if $h(x)$ gets small much faster than $g(x)$
out in the tails.  In such a case, though it is improbable (by definition) that
you will realize a value $x^{(i)}$ from the far tails of $h(x)$, if you do, then your
Monte Carlo estimate will be very large---$g(x^{(i)})/h(x^{(i)})$ for such an
improbable
$x^{(i)}$ may be orders of magnitude larger than the typical values $g(x)/h(x)$ that
you see.  On the other hand, if no or few such values $x^{(i)}$ are realized from
the far tails of $h(x)$, then the Monte Carlo estimate tends to be too small,
which is equally undesirable.   Such cases make importance sampling difficult, and
underscore the importance of choosing a good importance sampling function $h(x)$.


\subsection{Rao-Blackwellized estimators\label{sec:RBIntro}}
Especially in  the Bayesian analysis of complex stochastic models, one may be
interested in approximating the marginal posterior probability distributions of many
different quantities of interest.  In performing Monte Carlo in such models, many
different variables are simulated, and by using or combining those
simulated variables in different ways, one can derive different Monte Carlo
estimators for the same quantity.  Some of these estimators will be preferable to
others.  A variance reduction technique named ``Rao-Blackwellization"
\cite{Liuetal1994} provides a guideline for determing which of the possible Monte
Carlo estimators for a quantity should be used.
 I describe Rao-Blackwellization here in the context in which it is
employed throughout the dissertation---in the Monte Carlo estimation of
probabilities.  

Intuitively, the principle of Rao-Blackwellization can be understood as follows: 
let $W$ be a random variable describing the probability that a variable, say
$\lambda$, falls in a set
$\mathcal{L}$.  Then, assuming we could simulate values of $W$, we could estimate
$P(\lambda\in\mathcal{L})$ by the Monte Carlo estimate, $\frac{1}{n}\sum_{i=1}^n
w^{(i)}$.  Each $w^{(i)}$ will be a value between zero and one.  An alternative
estimate could be proposed: $\frac{1}{n}\sum_{i=1}^n
a^{(i)}$, where each $a^{(i)}$ takes the value 0 if
a uniform real number on $(0,1)$, say $U^{(i)}$ is less than $w^{(i)}$, and the
value 1 if $U^{(i)} > w^{(i)}$.  Doing so would be naive, because some information
is lost in condensing the real-valued $w^{(i)}$'s into the integer-valued
$a^{(i)}$'s.  Nonetheless, this sort of loss of information is routinely carried
out by people doing complicated Monte Carlo studies.  Rao-Blackwellization, in the
limited context I will use it in the dissertation, is the name given to improving
upon a Monte Carlo estimate of the form
$\frac{1}{n}\sum_{i=1}^n a^{(i)}$ by using, instead, a Monte Carlo estimate
of the form $\frac{1}{n}\sum_{i=1}^n
w^{(i)}$.  The general context in which it arises in the following chapters is
described below.  


Let us suppose that we have a probability model with three variables
(which, again, may be multidimensional), $X,Y$ and $\lambda$, having the joint
distribution $P(X,Y,\lambda)$.  Suppose that the value of $Y$ is known, and the
quantity that we desire to know is $P(\lambda\in\mathcal{L}|Y)$, the conditional
probability, given $Y$, that $\lambda$ is in some set $\mathcal{L}$.  This quantity
may be written as 
\begin{eqnarray}
P(\lambda\in\mathcal{L}|Y) &=& \sum_{x\in\mathcal{X}} P(\lambda\in\mathcal{L},X=x|Y)
\nonumber \\
& = & \sum_{x\in\mathcal{X}} P(\lambda\in\mathcal{L}|Y,X=x)P(X=x|Y)
\nonumber \\  &=&
\Exp[P(\lambda\in\mathcal{L}|Y,X)| Y].
\end{eqnarray}
Thus, a Monte Carlo estimator of $P(\lambda\in\mathcal{L}|Y)$ may
 be obtained as 
\begin{equation}
P(\lambda\in\mathcal{L}|Y) \approx
\frac{1}{n}\sum_{i=1}^n P(\lambda \in\mathcal{L}|Y,X^{(i)})
\label{eq:RB}
\end{equation}
where $X^{(i)}$ is simulated from $P(X|Y)$.  This would, in fact, be the
Rao-Blackwellized estimator for $P(\lambda|Y)$. Nonetheless, its use is sometimes
overlooked, for the following reason:  in many cases such as this, it is impossible
to obtain samples
$X^{(i)}$ ``directly" from $P(X|Y)$.  Rather, only samples
$(X^{(i)},\lambda^{(i)})$,
$i=1,\ldots,n$, are available, and it is thus tempting to apply (\ref{eq:ProbMC})
(rather than using (\ref{eq:RB})) to obtain the Monte Carlo estimator 
\begin{equation}
P(\lambda\in\mathcal{L}|Y) \approx
\frac{1}{n}\sum_{i=1}^n \ind\{\lambda^{(i)}\in\mathcal{L}\}.
\label{eq:notRB}
\end{equation}
This is often a poor choice, however, because, for independent samples
$(\lambda^{(i)},X^{(i)},Y^{(i)})$, the variance of (\ref{eq:notRB}) can be shown to
always be greater than or equal to the variance of the Monte Carlo estimator of
(\ref{eq:RB}) \cite{Gel&Smi1990}.  This may be proved using the Rao-Blackwell theorem
(hence the name) applied to the case where
$P(\lambda\in\mathcal{L}|Y)$ is the ``parameter" to be estimated, (\ref{eq:notRB})
is the unbiased estimator  based on an insufficient statistic for the parameter, and
(\ref{eq:RB}) is the Rao-Blackwellized version of (\ref{eq:notRB}).  

When the samples $(\lambda^{(i)},X^{(i)},Y^{(i)})$ are dependent across $i$, as is
the case with Markov chain Monte Carlo,  the superiority of the estimator
(\ref{eq:RB}) is more difficult to establish.  \citeN{Liuetal1994} prove the
superiority of (\ref{eq:RB}) for estimating expectations of functions of either
$\lambda$ or $X$, alone, with Markov chains having certain properties.  However,
they also present an example in which a Rao-Blackwellized estimator has higher
variance than the simple estimator.  Nonetheless, experience suggests that in most
cases involving the MCMC estimation of probabilities, the Rao-Blackwellized
estimator will perform better, especially as regards estimating probabilities in the
tails of the distribution.   

Quite often, the quantity  $P(\lambda \in\mathcal{L}|Y,X^{(i)})$ must be computed
anyway during the process of realizing a simulated pair of variables
$(X^{(i)},\lambda^{(i)})$.  In such cases the variance reduction of (\ref{eq:RB})
comes ``for free."  The practitioner of Monte Carlo should always be on the lookout
for opportunities to Rao-Blackwellize  Monte Carlo estimators. 


\subsection{Markov chain Monte Carlo---using dependent samples of $X$}
The ``Monte Carlo" part of Markov chain Monte Carlo is essentially identical to
regular Monte Carlo (except assessing convergence, \etc).  The main difference is
that the elements of the Monte Carlo sample $(X^{(1)},\ldots,X^{(n)})$ are not
independent of one another.  Rather they are sampled as states visited by an ergodic
Markov chain.  The main novelty is the need for some method to construct a Markov
chain with the appropriate limit distribution.  Almost all techniques for doing so
are variants of the Metropolis-Hastings sampler \cite{Hastings1970}.  This method
provides a way of constructing a time-reversible Markov chain with limit
distribution $\pi$  by satisfying, at each transition of the chain, the
detailed balance condition implied by $\pi$. This may be achieved when $\pi$ is
known only up to a normalizing constant.  Details of this may be found in numerous
texts.   





