# The Ball Pit Algorithm: A Markov Chain Monte Carlo Method Based on Path
Integrals

Miguel Fudolig
Department of Statistics, University of Nebraska-Lincoln
and
Reka Howard
Department of Statistics, University of Nebraska-Lincoln  Corresponding
Author<EMAIL_ADDRESS>
###### Abstract

The Ball Pit Algorithm (BPA) is a novel Markov chain Monte Carlo (MCMC)
algorithm for sampling marginal posterior distributions developed from the
path integral formulation of the Bayesian analysis for Markov chains. The BPA
yielded comparable results to the Hamiltonian Monte Carlo as implemented by
the adaptive No U-Turn Sampler (NUTS) in sampling posterior distributions for
simulated data from Bernoulli and Poisson likelihoods. One major advantage of
the BPA is its significantly lower computational time, which was measured to
be at least 95% faster than NUTS in analyzing single parameter models. The BPA
was also applied to a multi-parameter Cauchy model using real data of the
height differences of cross- and self-fertilized plants. The posterior medians
for the location parameter were consistent with other Bayesian sampling
methods. Additionally, the posterior median for the logarithm of the scale
parameter obtained from the BPA was close to the estimated posterior median
calculated using the Laplace normal approximation. The computational time of
the BPA implementation of the Cauchy analysis is 55% faster compared to that
for NUTS. Overall, we have found that the BPA is a highly efficient
alternative to the Hamiltonian Monte Carlo and other standard MCMC methods.

Keywords: Bayesian methods, Bayesian analysis, Path integral formulation,
Markov Processes, MCMC

##  1 Introduction

One of the key points in performing Bayesian analysis is to sample from the
posterior distribution of the model parameters. The posterior distribution,
which is the probability density function of the parameter, $\theta$,
conditional on the data $X$ can be calculated using the Bayes’ rule given by,

$P(\theta|X)=\frac{P(X|\theta)P(\theta)}{P(X)},$ (1)

where $P(X)$ is the marginal distribution of the data, $P(X|\theta)$ is the
likelihood of the data conditional on the parameter, and $P(\theta)$ is the
prior distribution of the parameter. The prior distribution P($\theta$)
incorporates any prior belief or knowledge about the parameter while the
likelihood accounts for the contribution of the data to the value of the
parameter (Robert, 2007). The posterior distribution is used to make
predictions, estimations, and decisions in the Bayesian framework.

There are special cases where the posterior distribution is described by a
standard, well-defined distribution and are thus easy to be sampled from. A
good example of these distributions are conjugate priors, which is a group of
distributions where the resulting posterior distribution is in the same family
of distributions as the defined prior (Gelman et al., 2013). One example of a
conjugate prior is the beta prior for a binomial likelihood where evaluating
Equation 1 when the prior $P(\theta)$ is a beta distribution and the
likelihood function $P(X|\theta)$ is a binomial distribution yields a beta
distribution for the posterior distribution. An advantage of using conjugate
priors is the ease of sampling from standard distributions.

However, not all likelihoods have a conjugate prior distribution and hence, we
need to resort to numerical methods for sampling posterior distributions such
as Markov chain Monte Carlo (MCMC) methods. MCMC methods use Markov chains of
sequential draws of the parameter from approximate distributions to get a good
approximation of the target posterior distribution (Gelman et al., 2013). Some
well-known MCMC algorithms include the Gibbs sampler (Casella and George,
1992; Geman and Geman, 1984; Smith and Roberts, 1993), Sequential Monte Carlo
(SMC) method (Liu and Chen, 1998; Doucet et al., 2001; Del Moral, 1996), and
the Particle Markov Chain Monte Carlo (p-MCMC) method (Andrieu et al., 2009,
2010; Holenstein, 2009). The state-space models used in the aforementioned
MCMC algorithms simulate a Markov chain in the state space and introduce a
criterion for accepting proposal points from the Markov chain based on the
target distribution (Gelman et al., 2013; Yuan et al., 2012). The Hamiltonian
Monte Carlo (HMC) method (Betancourt, 2017; Duane et al., 1987; Hoffman et
al., 2014) improves the efficiency of sampling the target distributions by
combining MCMC with deterministic simulation methods through incorporating the
concept of the Hamiltonian function from physics (Gelman et al., 2013).
Girolami and Calderhead (2011) improved the efficiency of the HMC by proposing
the Riemann Manifold Hamiltonian Monte Carlo (RMHMC), which uses the Riemann
geometry of high-dimensional state-spaces to improve the sampling efficiency
of the HMC. Lan et al. (2015) decreased the numerical cost of the RMHMC by
proposing an implicit geometric integrator, which resulted to a Lagrangian
formulation of the RMHMC.

Another way to decide whether to accept a proposal point is to determine how
likely will the chain continue on the trajectory defined by the proposal
point. This information can be provided by path integrals, which employs the
concept that the total transition probability between two states A and B is
the sum of the probabilities of all possible paths from state A to state B
(Feynman et al., 2010; Wio, 2013; Albeverio et al., 1976). Path integrals are
widely used in quantum and statistical physics in investigating time evolution
of physical systems (Feynman et al., 2010; MacKenzie, 2000; Wio, 2013;
Hawking, 1979; Caldeira and Leggett, 1983; Wiegel, 1975), but the concept of
path integrals have not been widely utilized in Bayesian analysis. Path
integrals were used by Chang et al. (2014) to perform Bayesian inference on
physical inverse problems. Chang et al. (2015) used path integrals to develop
a non-parametric Bayesian approach to quantify bond energies and mobilities.
However, these aforementioned papers did not apply the path integral
formulation to Markov chains in state-space models. Fujii and Hatakenaka
(2018) applied the path integral approach to derive analytic expressions of
posterior distributions of Markov chains with a specific form of likelihood.
However, the derived expressions were not utilized in estimating or making
inference on these parameters.

Applying the concept of path integrals to Markov chains in state-space models
gives us crucial information about the path taken by the Markov process, which
can be used to devise a faster algorithm in sampling posterior distributions.
The resulting algorithm is expected to be more efficient compared to the
common MCMC methods because the solution of the path integral provides
information that predicts the trajectory of the Markov chain in the state-
space, which we also refer to as the parameter space, without the need to
sample momentum values at each time step.

We developed a new Markov Chain Monte Carlo algorithm to approximate marginal
posterior distributions using the resulting equation from the path integral
formulation for a Wiener process, a special type of Markov process that
describes Brownian motion (Szabados, 2010), in the parameter space. This
algorithm was applied on simulated data for Bernoulli and Poisson models and
the estimated summary statistics of the posterior distribution was compared to
the estimates obtained by the R package rstan (R Core Team, 2017; Stan
Development Team, 2020). The package rstan implements the Hamiltonian Monte
Carlo using an adaptive approach of the No U-Turn Sampler (NUTS) (Hoffman et
al., 2014). The resulting algorithm was then applied to a multiparameter
Cauchy model with real-world data from the height differences of cross- and
self-fertilized plants from the R package LearnBayes (R Core Team, 2017;
Albert, 2018).

Section 2 provides the theoretical basis of the algorithm by introducing the
concept of path integrals and explaining its application to the Bayesian
analysis of Markov chains. Section 3 explains the proposed algorithm, which we
call the Ball Pit Algorithm (BPA), and show how it compares against results
from NUTS. This is followed by a discussion of the results after analyzing the
simulated data using the BPA and comparing it to the NUTS results.

##  2 The Path Integral Formulation

The path integral method has been used in understanding the evolution of
physical systems in time. In quantum mechanics, the path integral can be used
in determining how the state of a particle evolves in time by accounting its
interactions with its environment (Feynman et al., 2010; Sakurai and
Napolitano, 2011). Wio (2013) worked on the application of path integrals to
stochastic processes such as Wiener processes. In this section, we apply the
concept of path integrals to Bayesian analysis after briefly introducing what
path integrals are.

###  2.1 Path Integral

Consider points A at position $\mathbf{x}$ and B at position $\mathbf{y}$. If
one particle starts at point A at time $t_{0}$ and ends up at point B at time
$t$, there are various paths that the particle could have taken from point A
to reach point B in between time $t_{0}$ and $t$. Some of these paths are
shown in Figure 1. The path integral formulation assumes that the total
transition probability $P(\mathbf{y},t|\mathbf{x},t_{0})$ is the sum of the
probabilities of all possible paths the particle could have taken from point A
to point B.

Figure 1: The blue, green, and orange paths are some of the paths that a
particle could take from point A to point B from time $t_{0}$ to $t$.

According to Feynman et al. (2010), the conditional probability
$P(\mathbf{y},t|\mathbf{x},t_{0})$ that the system can have a value of
$\mathbf{y}$ at time $t$ given that it had a value of $\mathbf{x}$ at time
$t_{0}$ can be solved by defining a sequence of values $\\{x_{k}\\}_{k=1}^{N}$
such that $x_{0}=x$, $x_{N}=y$, and $t=t_{0}+N\epsilon$ where $\epsilon>0$ is
the time increment. The resulting formula of the conditional probability is:

$\displaystyle P(\mathbf{y},t|\mathbf{x_{0}},t_{0})$ $\displaystyle=\idotsint
dx_{1}dx_{2}...dx_{N-1}P(y,t^{\prime}|x_{N-1},t+\epsilon)$
$\displaystyle\cdot\cdot\cdot
P(x_{2},t+2\epsilon|P(x_{1},t+\epsilon)P(x_{1},t+\epsilon|x_{0},t)P(x_{0},t).$
(2)

This can then be written in terms of the path integral notation as,

$P(\mathbf{y},t|\mathbf{x},t_{0})=\int\mathcal{D}x[t]\Phi[x(t)]$ (3)

where

$\mathcal{D}x[t]=\lim_{n\to\infty}dx_{1}dx_{2}...dx_{n-1}$ (4)

and

$\Phi[x(t)]=\lim_{\epsilon\to
0}\prod_{i=1}^{N}P(x_{i},t+i\epsilon|x_{i-1},t+(i-1)\epsilon)$ (5)

The term $\Phi[x(t)]$ is referred to as the propagator (Sakurai and
Napolitano, 2011), which can be interpreted as a probability of a specific
path the system follows from $x_{0}$ to $x_{N}$. For stochastic processes such
as Brownian motion and Markov chains, Equation 3 can be simplified into the
following form,

$P(\mathbf{y},t|\mathbf{x},t_{0})=\int\mathcal{D}x[t]\exp[-S_{E}]$ (6)

where $S_{E}$ is the discrete-time Euclidean action can be defined as the
following,

$S_{E}=\sum_{i=1}^{N}\mathcal{L}(x_{i},\dot{x}_{i},t_{i})\Delta t_{i}$ (7)

where the integration runs over time from $t=t_{1}$ to $t^{\prime}=t_{N}$ and
$\dot{x}=\dfrac{x(t)-x(t-\epsilon)}{\epsilon}$ is the speed of the system. The
quantity $\mathcal{L}$ is defined as the Lagrangian of the system, which is a
key function in physics in describing the motion of the system (Goldstein et
al., 2002). The Lagrangian of the system is expressed as the difference
between the kinetic energy $K$ and the potential energy $U$ (Goldstein et al.,
2002). The kinetic energy is associated with the energy in motion while the
potential energy is associated with the configuration of the system (Tipler
and Mosca, 2007). Once the Lagrangian function of a process is established,
the trajectory of the process is well-defined given that its initial position
and velocity is known.

We can apply these concepts to Bayesian analysis by utilizing path integrals
to describe the trajectory of Markov chains used in MCMC methods. The path
integral formulation of Bayesian analysis of Markov chains is explained
further in the next subsection.

###  2.2 Path Integrals: Bayesian Analysis of Markov Chains

Fujii and Hatakenaka (2018) used the path integral approach for Markov chains
in Bayesian analysis. One of the limitations of their study was that their
approach only applied to likelihoods with the form:

$L(\theta|x)=M\exp\left[-V(\theta-x(t))\epsilon\right],$ (8)

where $\theta$ is the parameter, $x(t)$ is the data, $\epsilon$ is the time
increment between the steps of the Markov chain, $V$ is a function of the
difference of the parameter and the data, and $M$ is the normalization
constant. This form of the likelihood is similar to exponential family
likelihood except for the explicit dependence on the time increment $\epsilon$
and the difference between parameter and data, which are not commonly
encountered in statistical models. We derive the path integral formulation for
a more general form of the likelihood and study the implications of the
resulting expression.

We start with the concept of recursive Bayesian estimation for Markov chains
outlined by Chen et al. (2003). Suppose we have a parameter $\theta_{N}$ at
time step $N$ that follows a first-order Markov process, i.e. $\theta_{N}$ is
only conditional on the value of the immediately preceding time step and
independent of the rest of the steps in the path. We denote
$\mathbf{X}_{N}=(x_{0},x_{1},...,x_{N})$ as the vector of independent
observations from time-step $t=0$ to $t=N$. Then the posterior distribution of
$\theta_{N}$ conditional on the observations $\mathbf{X}_{N}$ can be solved
using the Bayes’ rule.

$\displaystyle P(\theta_{N}|\mathbf{X}_{N})$ $\displaystyle=$
$\displaystyle\dfrac{P(\mathbf{X}_{N}|\theta_{N})P(\theta_{N})}{P(\mathbf{X}_{N})}$
(9) $\displaystyle=$
$\displaystyle\dfrac{P(x_{N},\mathbf{X}_{N-1}|\theta_{N})P(\theta_{N})}{P(x_{N},\mathbf{X}_{N-1})}$
(10) $\displaystyle=$
$\displaystyle\dfrac{P(x_{N}|\theta_{N},\mathbf{X}_{N-1})P(\mathbf{X}_{N-1}|\theta_{N})P(\theta_{N})}{P(x_{N}|\mathbf{X}_{N-1})P(\mathbf{X}_{N-1})}$
(11) $\displaystyle=$
$\displaystyle\dfrac{P(x_{N}|\theta_{N})P(\theta_{N}|\mathbf{X}_{N-1})}{P(x_{N}|\mathbf{X}_{N-1})}$
(12)

The term $P(\theta_{N}|\mathbf{X}_{N-1})$ is the probability density of the
parameter conditioned on the observations before the current time step, which
implies that this term is a conditional prior for $\theta_{N}$. This
probability can be expressed as the following integral:

$P(\theta_{N}|\mathbf{X}_{N-1})=\int
P(\theta_{N}|\theta_{N-1})P(\theta_{N-1}|X_{N-1})d\theta_{N-1},$ (13)

Substituting 13 into 12,

$P(\theta_{N}|\mathbf{X}_{N})=\dfrac{P(x_{N}|\theta_{N})}{P(x_{N}|\mathbf{X}_{N-1})}\int
P(\theta_{N}|\theta_{N-1})P(\theta_{N-1}|\mathbf{X}_{N-1})d\theta_{N-1}$ (14)

where the integral is on the support of $\theta_{N-1}$. The term
$P(\theta_{N}|\theta_{N-1})$ can be interpreted as the transition probability
of the parameter from one time step to another. We consider the parameter
$\theta$ as a random variable that evolves as described by the Brownian motion
stochastic differential equation given by:

$d\theta=\mu(\theta,t)dt+\sigma dW_{t},$ (15)

where $\mu(\theta,t)$ is the drift velocity of the stochastic process and
$\sigma$ is the strength of the noise. $W_{t}$ is known as the Wiener process
at time $t$, which is the continuous model of Brownian motion. The Wiener
process can also be defined in terms of the following distribution (Szabados,
2010; Wio, 2013),

$W_{t}-W_{s}\sim N(0,s),$ (16)

where $s$ is a time point such that $0\leq s<t$. When the drift velocity is
set to zero, Equation 15 reduced to a Wiener process for $\theta$. Based on
the definition of the Wiener process, the distribution of $\theta_{t+dt}$ for
any time change $dt>0$,

$\theta_{t+dt}|\theta_{t}\sim N(\theta_{t},\sigma^{2}dt)$ (17)

This can also be written in equation form,

$P(\theta_{t+dt}|\theta_{t})=(2\pi\sigma^{2}dt)^{-1/2}e^{-\frac{\dot{\theta}_{t+dt}^{2}}{2\sigma^{2}}dt}$
(18)

where $\displaystyle\dot{\theta}_{t+dt}=\frac{\theta_{t+dt}-\theta_{t}}{dt}$.

We consider $dt=1$, which corresponds to the discrete case of the Markov
process. Substituting Equation 18 to Equation 14 yields the following
integral:

$P(\theta_{N}|\mathbf{X}_{N})=\dfrac{P(x_{N}|\theta_{N})}{P(x_{N}|\mathbf{X}_{N-1})(2\pi\sigma^{2})^{1/2}}\int{e^{-\dfrac{(\theta_{t+1}-\theta_{t})^{2}}{2\sigma^{2}}}P(\theta_{N-1}|\mathbf{X}_{N-1})d\theta_{N-1}}$
(19)

This can also be written as,

$P(\theta_{N}|\mathbf{X}_{N})=C\int{e^{\left[-\dfrac{(\theta_{N}-\theta_{N-1})^{2}}{2\sigma^{2}}+\log
P(x_{N}|\theta_{N})\right]}P(\theta_{N-1}|\mathbf{X}_{N-1})d\theta_{N-1}},$
(20)

where all the terms constant with respect to $\theta_{N}$ are pooled into a
general normalizing constant $C$. Equation 20 shows a recursive pattern
between the posterior distributions at time steps $N-1$ and $N$. Continuing
this recursion, we get the following integral,

$P(\theta_{N}|\mathbf{X}_{N},\theta_{0})=C\idotsint\prod_{i=0}^{N-1}d\theta_{i}e^{\left[-\dfrac{(\theta_{i+1}-\theta_{i})^{2}}{2\sigma^{2}}+\log
P(x_{i+1}|\theta_{i+1})\right]}P(\theta_{0}),$ (21)

We focus on the product in the integrand which can be written as,

$\prod_{i=0}^{N-1}d\theta_{i}e^{\left[-\dfrac{(\theta_{i+1}-\theta_{i})^{2}}{2\sigma^{2}}+\log
P(x_{i+1}|\theta_{i+1})\right]}=\exp\left[-\sum_{i=1}^{N}\left(\dfrac{\dot{\theta}_{i}^{2}}{2\sigma^{2}}-\log
P(x_{i}|\theta_{i})\right)\right]$ (22)

The above equation is similar to the definition of the Wiener path integral
related to the Wiener measure (Glimm and Jaffe, 2012; Wio, 2013). Comparing
the terms to Equation 6, we determine the discrete Euclidean action for the
parameter space of $\theta$ is given by,

$S_{E}=\sum_{i=1}^{N}\left(\dfrac{\dot{\theta}_{i}^{2}}{2\sigma^{2}}-\log
P(x_{i}|\theta_{i})\right)$ (23)

Based on Equation 23 and the idea that the process should accelerate towards
the point of maximum likelihood, the corresponding Lagrangian to determine the
equation of motion of the process is given by,

$\mathcal{L}=\left[\dfrac{\dot{\theta}_{\tau}^{2}}{2\sigma^{2}}+\left(\log
P(x_{\tau}|\theta_{\tau})\right)\right]$ (24)

Classical physics tells us that the process will evolve by taking the path of
least action. Since the Euclidean action is dependent on the likelihood, this
implies that points simulated from the prior will evolve in the parameter
space as dictated by the value of the log-likelihood. As we reach equilibrium,
the prior distribution is expected to have evolved to the posterior
distribution as guided by the likelihood.

The first term of Equation 24 can be interpreted as the energy due to the
inherent motion of the process due to its stochastic nature. A process with a
high $\sigma$ value will have a larger variance in its trajectory and hence
will cover a wider area in the parameter space. This is analogous to the
smaller masses accelerating faster when a force is applied to them, which is
consistent with the interpretation of noise strength by Lan et al. (2015) in
their work on Lagrangian Monte Carlo.

The log-likelihood in Equation 23 incorporates the information from the
observations on choosing the path taken by the process. We expect the process
to accelerate towards the value of maximum likelihood or equivalently, the
value of the minimum negative log-likelihood. This implies that the process
should be pulled towards this point at any given time. Using the negative log-
likelihood as the counterpart for the potential energy instead of the positive
makes the dynamics of the process consistent with what is expected in physical
systems. We use this idea in formulating the equations of motion in the
parameter space.

###  2.3 Equations of Motion

The Lagrangian can be used to find the equations of motion of the system
(Goldstein et al., 2002) which, in the context of this research, will describe
how the system moves in the parameter space upon recording an observation,
$x$, at a given time $\tau$. Assuming that our likelihood is explicitly
independent of the speed of the process, $\dot{\theta}$, and $\tau$, the
equations of motion for the parameter is given by the Euler-Lagrange equation
(Goldstein et al., 2002),

$\dfrac{d}{dt}\left(\dfrac{\partial\mathcal{L}}{\partial\dot{\theta}}\right)-\dfrac{\partial\mathcal{L}}{\partial\theta}=0$
(25)

This principle is similar to the Hamiltonian equations of motion that we use
in Hamiltonian Monte Carlo (Betancourt, 2017), but instead of sampling
momentum values for every time step, the Euler-Lagrange equation can be used
to update the speed of the process in the state space. This result describes
how the process moves in the $\theta$-space, which can be interpreted as the
parameter space in the context of Bayesian analysis. The form of Equation 24
is similar to the Lagrangian function used by Lan et al. (2015) in their work
on Explicit Riemannian Manifold Lagrangian Monte Carlo. The resulting equation
is a second order ordinary differential equation (ODE) in time which is easily
solvable using numerical methods.

We use this form for the equation of motion of the process to create an
algorithm that samples from a target posterior distribution by evolving points
sampled from the prior in an MCMC framework.

##  3 The Ball Pit Algorithm (BPA): An Approximate Bayesian Method to Estimate
Posterior Statistics

In this section, we discuss how a Monte Carlo algorithm to approximate
posterior distributions was developed using the Lagrangian function and the
equation of motion we derived in Section 2.2.

###  3.1 Algorithm

The main idea of the algorithm is that the prior distribution evolves to the
posterior distribution by taking the path dictated by the value of the
likelihood. A Wiener process moving in the parameter space will follow the
trajectory dictated by the equation of motion derived from the Lagrangian
function.

A pseudocode of the algorithm is presented in Algorithm 1. The algorithm
starts by simulating starting points from the prior distribution, an approach
used by Pompe et al. (2020) in their adaptive MCMC algorithm for multimodal
distributions. These points will be utilized as starting points for a single
independent process which we will refer to as a ”ball”. Each ball will be
given an initial speed, $\dot{\theta}$, which will be simulated from a
Gaussian distribution with centered at zero and variance equal to the noise
strength of the Markov process $\sigma^{2}$. The trajectory of each ball,
which we refer to as the chain, in the state space will be recorded by
numerically solving the differential equation in Equation 25. In the evolution
of these chains, we need to avoid moving on paths unlikely for the chain to
follow. Hence, we need to introduce a stopping criterion such that when the
criterion is satisfied, the ball would change its trajectory to avoid
following an unlikely path. This stopping criterion will be based on the
probability of following the proposed path, $E$, given by,

$E(\theta,\dot{\theta})=\exp\left[{-\epsilon\frac{\dot{\theta}^{2}}{2\sigma^{2}}+\log
L(X|\theta)}\right],$ (26)

where $\epsilon$ is the time increment between steps in solving the equation
of motion and $\log L(X|\theta)$ is the log-likelihood. In the context of the
path integral formulation, the effective weight of the proposed path is the
propagator. We impose the stopping criterion for a candidate point
$(\theta^{*},\dot{\theta}^{*})$ given a current point
$(\theta_{0},\dot{\theta}_{0})$ in a form of an acceptance probability,
$p_{accept}$, defined as:

$p_{accept}=\min\left(1,\frac{E(\theta^{*},\dot{\theta^{*}})}{E(\theta_{0},\dot{\theta_{0}})}\right)$
(27)

The term
$\displaystyle\frac{E(\theta^{*},\dot{\theta^{*}})}{E(\theta_{0},\dot{\theta_{0}})}$
can be interpreted as a risk ratio in continuing the movement along the
proposed path. The candidate point will be accepted with probability
$p_{accept}$. If the candidate point is rejected, the process will remain at
$\theta_{0}$ and its speed will be re-sampled from the same Gaussian
distribution where the initial speeds were sampled from.

There are situations where the value for $\theta^{*}$ is rejected for each
step depending on the choice for the prior distribution. To avoid the ball
being stuck at an unlikely parameter value, we perform a check for every
rejection if the point has been rejected for a user-defined lag time $\tau$.
If this is the case, the ball will be given a new starting position by
sampling from the prior distribution again.

All balls are allowed to move until a final evolution time $T$ when the
desired chain length is reached. After accounting for a burn-in phase for each
chain, the chains are combined to form an approximate posterior distribution.
We call this algorithm the Ball Pit Algorithm (BPA) because the posterior
distribution is approximated by the distribution of the balls, which are
assumed to cover the relevant area of the parameter space.

This algorithm was tested on single parameter models and the posterior
statistics obtained from the BPA will be compared to the results of the
adaptive No U-Turn Sampler (NUTS) algorithm, which is commonly used to
implement HMC. The BPA was also used to analyze real-world data of plant
heights of self- and cross-fertilized plants by sampling marginal posterior
distributions for a Cauchy multi-parameter model.

Algorithm 1 The Ball Pit Algorithm: One Process

1:procedure Path Integrals

2: $\text{Initialize}~{}\textit{prior,}~{}\sigma^{2},\epsilon,t=0,T,\tau$

3: $\theta_{0}\sim\textit{prior}$

4: $\dot{\theta_{0}}\sim N(0,\sigma^{2})$

5:  _Chain Start_ :

6: $\theta\leftarrow\theta_{0}$

7: $v\leftarrow\dot{\theta_{0}}$

8:  _Update_ :

9: $\text{Update}~{}(\theta,v)~{}\text{using Euler-Lagrange Equations to
get}~{}(\theta^{*},v^{*})$

10: if $L(\theta^{*})~{}\text{is NA}$ then

11: $\text{keep }\theta$

12: $v\sim N(0,\sigma^{2})$

13:
$p_{acc}\leftarrow\min(1,\exp{\left[E(\theta^{*},v^{*})-E(\theta,v)\right]})$

14: $p\sim Unif(0,1)$

15: if $p<p_{acc}$ then

16: $\theta\leftarrow\theta^{*}$

17: $v\leftarrow v^{*}$

18: else

19: if $\theta=\theta_{(t-\tau)}$ then

20: $\text{Resample }\theta\text{ from prior.}$

21: $v\sim N(0,\sigma^{2})$

22: else

23: $\text{keep}~{}\theta$

24: $v\sim N(0,\sigma^{2})$

25: $t\leftarrow t+\epsilon$

26: if t = T then

27: End Process

28: else

29: go to _Update_.

###  3.2 Advantages

The proposed method is similar to the Hamiltonian Monte Carlo, which is also a
method derived from physics concepts. The Hamiltonian Monte Carlo (HMC)
involves simulating ”momentum” values, $p$, at the start of each iteration
from a conditional normal distribution to determine the gradient of steepest
descent of the target distribution in the phase space ($(\theta,p)$ space).
The proposed value of the parameter is obtained from the Hamiltonian equation
of motion (Betancourt, 2017). The integration of physics in the HMC decreases
the noise due to the inherent stochasticity of the Markov chain, making it
more efficient compared to the other MCMC methods (Gelman et al., 2013).

In contrast, the BPA follows a Lagrangian framework where the position and the
velocity of the process evolves as dictated by the Euler-Lagrange equation
(Equation 25). New velocities for each ball need to be sampled only when the
chain starts or the proposed point is rejected, which means there are fewer
sampling steps in the BPA compared to the HMC. This translates to faster
convergence of the algorithm and lower computational time compared to the HMC.

##  4 Bayesian Analysis Using the BPA: Simulated and Real-World Data

This section discusses the results obtained from using the Ball Pit Algorithm
(BPA) to approximate the posterior distribution of single parameter models. In
this research, we considered single parameter models with Bernoulli and
Poisson likelihoods. The Bernoulli likelihood was chosen because of its
parameter has a bounded support and the existence of a conjugate prior. The
Jeffreys prior for the Poisson distribution is an improper prior, which
enabled us to explore how improper priors can affect the resulting
approximation for the posterior distributions. Posterior statistics and the
computation time for the algorithm were recorded and then compared to the
results obtained through the NUTS algorithm included in the R package rstan
(Stan Development Team, 2020; R Core Team, 2017).

The BPA was implemented using 80 balls that are made to move around the state
space for an evolution time of 1 second with a time step size of 0.01 seconds
to match the default overall chain length used by rstan (Stan Development
Team, 2020). The values for the first five seconds, which translates to a
chain length of 500 for each ball, will be treated as a warm-up phase and will
be discarded just as in the regular MCMC algorithm. The sampled marginal
posterior distribution was then compared to a NUTS implementation that used 80
chains with a chain length of 100.

We then tested the BPA to sample marginal posterior distributions from a
multi-parameter model using a real-world data set from the LearnBayes package
used by Albert (2009) that contains the height differences between cross- and
self-fertilized plants.

###  4.1 Bernoulli Process

The corresponding likelihood for a sample size $N$ for a Bernoulli process
with probability $\theta$ can be written as,

$P(X|\theta)=\prod_{i=1}^{N}\theta^{x_{i}}(1-\theta)^{N-x_{i}}$ (28)

We used Equation 24 to derive the Lagrangian function with the given
likelihood function with probability $\theta$, size $N$ and observation
$\mathbf{x}=(x_{1},x_{2},...x_{N})$.

$\mathcal{L}=\left[\frac{\dot{\theta}^{2}}{2\sigma^{2}}+\left(\sum_{i=1}^{N}x_{i}\log\theta+\left(N-\sum_{i=1}^{N}x_{i}\right)\log(1-\theta)\right)\right]$
(29)

Then the corresponding equation of motion for $\theta$ can be solved from
Equation 25:

$\displaystyle\dfrac{d}{dt}\left(\dfrac{\dot{\theta}}{\sigma^{2}}\right)-\left[\dfrac{\sum_{i=1}^{N}x_{i}}{\theta}-\dfrac{\left(N-\sum_{i=1}^{N}x_{i}\right)}{1-\theta}\right]$
$\displaystyle=$ $\displaystyle 0$ (30)
$\displaystyle\ddot{\theta}-\sigma^{2}\left[\dfrac{\sum_{i=1}^{N}x_{i}-N\theta}{\theta(1-\theta)}\right]$
$\displaystyle=$ $\displaystyle 0$ (31)

where $\theta\in(0,1)$. $\dot{\theta}$ and $\ddot{\theta}$ are the velocity
and acceleration of the process, respectively. Solving the second-order
differential equation in Equation 31 would provide the candidate values for
the parameter $\theta$ and velocity $\dot{\theta}$ of the process.

For the simulation study, a representative data set of size 200 was simulated
from a Bernoulli distribution with probability set to a representative value
$\theta=0.3$. We set the Markov process variance to be 1. Two types of priors
will be considered for this analysis to see whether the consistency of the BPA
compared to the NUTS algorithm depends on the choice of prior. Figures 2 and 3
and Tables 2 show the respective approximate posterior distributions and
quantiles for the non-informative and informative prior for $\theta$. The non-
informative prior used is $Beta(1,1)$, which is the same as the $Unif(0,1)$
distribution, while the informative prior used is $Beta(3,7)$. The beta
distribution was selected to be the choice of prior because it is a conjugate
prior to a Bernoulli likelihood.

Figure 2: The samples from the posterior distribution from the BPA and NUTS
method for the uniform prior. The total number of post warm-up samples for
both methods is 4000. Figure 3: The samples from the posterior distribution
from the BPA and NUTS method for an informative Beta(3,7) prior. The total
number of post warm-up samples for both methods is 4000. Table 1: Estimates of
the posterior distribution from theoretical Bayes, NUTS, and the Ball Pit
Algorithm (BPA) with time elapsed for a uniform prior ($Beta(1,1)$).

Method Mean Std. Error. Std. Dev. 2.5 % 25 % 50 % 75 % 97.5 % Time (s) NUTS
0.30 <0.01 0.03 0.24 0.28 0.3 0.32 0.37 34.23 BPA 0.30 <0.01 0.03 0.25 0.28
0.3 0.32 0.36 0.71

Table 2: Estimates of the posterior distribution from NUTS and the Ball Pit
Algorithm (BPA) with time elapsed for an informative prior Beta(3,7).

Method Mean Std. Error. Std. Dev. 2.5 % 25 % 50 % 75 % 97.5 % Time NUTS 0.3
<0.01 0.03 0.24 0.28 0.30 0.32 0.37 29.39 BPA 0.3 <0.01 0.03 0.25 0.28 0.3
0.32 0.37 0.63

The two methods yielded similar approximations of the posterior distribution
as displayed in the plots in Figures 2 and 3. The quantile estimates of the
BPA are also close to the results provided by NUTS. The big difference between
the two methods lie in the computational time. The computational time was
reduced by 98% in generating the approximate posterior distribution for the
uniform prior, while the computational time was reduced by 97% for the
informative prior.

###  4.2 Poisson Process

If the response follows a Poisson distribution with a constant rate parameter
$\lambda$, then the likelihood can be written as

$P(X|\theta)=\prod_{i=1}^{N}e^{-\lambda}\frac{\lambda^{x_{i}}}{x_{i}!}.$ (32)

The corresponding Lagrangian is given by,

$\mathcal{L}=\left[\dfrac{\dot{\lambda}^{2}}{2\sigma^{2}}+\left(-\lambda+x\log\lambda-\log\left(x!\right)\right)\right]$
(33)

The equation of motion based on this Lagrangian can be calculated by taking
the logarithm of this likelihood and substituting to Equation 24,

$\displaystyle\dfrac{d}{dt}\left(\dfrac{\dot{\lambda}}{\sigma^{2}}\right)-\left[\dfrac{x}{\lambda}-1\right]$
$\displaystyle=$ $\displaystyle 0$ (34)
$\displaystyle\ddot{\lambda}-\sigma^{2}\left[\dfrac{x-\lambda}{\lambda}\right]$
$\displaystyle=$ $\displaystyle 0$ (35)

Similar to Equation 31, the solution of Equation 35 is well-defined if the
initial values of the rate parameter value $\lambda$ and velocity
$\dot{\lambda}$ are set.

For the simulation study, a data set of size 200 was simulated from a Poisson
distribution with a rate parameter $\lambda=40$. We used Equation 35 to update
the position and the speed of the process in the BPA. The value of the Markov
process variance was set to $\sigma^{2}=100$. We chose a larger value for
$\sigma^{2}$ because we want the process to have less ”mass” to cover a wider
area of the parameter space for the rate parameter, which is larger compared
to the parameter space of the Bernoulli probability in the previous example.
Similar to the approach we used for the Bernoulli likelihood example, we now
compare the performance of the BPA compared with NUTS for non-informative and
informative priors on a simulated Poisson data set. The Jeffreys prior was
used as the non-informative prior for the Poisson rate parameter, $\lambda$,
which is given by $P(\lambda)\propto\lambda^{-1/2}$. Meanwhile, the
informative prior used was a normal distribution centered at the sample mean
with a variance of 4. The density plots of the approximate distributions for
both priors are shown in Figures 4 and 5, and the quantile estimates are
listed in Tables 3 and 4.

Table 3: Estimates of the posterior distribution from theoretical Bayes, NUTS,
and the Ball Pit Algorithm (BPA) with time elapsed for the Jeffreys prior
$P(\lambda)\propto\lambda^{-1/2}$.

Method Mean Std. Error. Std. Dev. 2.5 % 25 % 50 % 75 % 97.5 % Time NUTS 40.35
0.01 0.46 39.45 40.04 40.36 40.64 41.26 37.41 BPA 40.34 0.08 0.40 39.58 40.04
40.34 40.64 41.08 1.17

Figure 4: The samples from the posterior distribution from the BPA and NUTS
method for the Jeffreys prior. The total number of post warm-up samples for
both methods is 4000. Table 4: Estimates of the posterior distribution from
NUTS and the Ball Pit Algorithm (BPA) with time elapsed for an informative
$N(\bar{x},4)$ prior.

Method Mean Std. Error. Std. Dev. 2.5 % 25 % 50 % 75 % 97.5 % Time NUTS 40.36
0.01 0.44 39.53 40.07 40.35 40.65 41.23 33.22 BPA 40.32 0.09 0.39 39.60 40.04
40.33 40.61 41.05 1.06

Figure 5: The samples from the posterior distribution from the BPA and NUTS
method for an informative $N(\bar{x},4)$ prior. The total number of post warm-
up samples for both methods is 4000.

The density plots for the posterior distribution show that the approximated
posterior distribution obtained using the BPA is narrower compared to the NUTS
result. This observation was supported by the smaller standard deviation and
narrower credible intervals for the BPA chain in Tables 3 and 4. The
difference in width between the two approximations was more pronounced for the
informative prior, which might be explained by the mechanism of the algorithm.
Since the starting points are sampled from the prior, the balls are localized
in the neighborhood of the minimum point of the negative log-likelihood. This
made it unlikely for the chains to explore points far from the mean, resulting
to a narrower distribution from the sampled points. Similar to the Bernoulli
likelihood example, the computational time was significantly reduced by 97%
for both priors.

###  4.3 Multiparameter Model: Cauchy Distributed Errors

We tested the BPA to sample marginal posterior distributions from a multi-
parameter model. We considered a real-world data set from the LearnBayes
(Albert, 2018) package used by Albert (2009) that contains the height
differences between cross- and self-fertilized plants. These differences are
modeled by a Cauchy likelihood with a location parameter $\mu$ and scale
parameter $\sigma$.

The approach we took in applying the BPA to a multiparameter model is similar
to Gibbs sampling where we make an educated guess about the parameter values
and sequentially updating these values through implementing the algorithm
recursively. We used the laplace function in the LearnBayes (Albert, 2018)
package to get the an estimate of the posterior modes for both the location
parameter and the logarithm of the scale parameter. We then use these values
as an input to the BPA as it simulates the marginal posterior distribution for
both parameters.

For the multiparameter model, the second-order partial differential equations
from the Euler-Lagrange equations were set up and solved numerically. The
joint prior used for this analysis is the non-informative prior $\displaystyle
P(\mu,\sigma)\propto\sigma^{-1}$. Both algorithms used 100 balls/chains with a
chain length of 1000 for their respective implementations. The resulting
approximate posterior distributions from both methods are shown in Figures 6
and 7.

Figure 6: The samples from the posterior distribution from the BPA and NUTS
method for the non-informative prior $\displaystyle
P(\mu,\sigma)\propto\sigma^{-1}$. The total number of post warm-up samples for
both methods is 50000. Figure 7: The samples from the posterior distribution
from the BPA and NUTS method for the non-informative prior $\displaystyle
P(\mu,\sigma)\propto\sigma^{-1}$. The total number of post warm-up samples for
both methods is 50000.

We see that there was a good correspondence between the shape of the posterior
distributions for both $\mu$ and $\log\sigma$, with the BPA yielding slightly
narrower results compared to the NUTS algorithm. Table 5 also shows the
resulting summary statistics for the two parameters from the two methods.

Table 5: Estimates of the posterior distribution from NUTS and the Ball Pit
Algorithm (BPA) for the multi-parameter Cauchy likelihood model.

Parameter Method Mean Std. Error. Std. Dev. 2.5% 25% 50% 75% 97.5% $\mu$ STAN
25.54 0.04 7.05 12.58 20.96 25.12 29.67 40.96 BPA 25.1 0.45 5.20 15.50 21.34
25.02 28.71 35.44 $\log\sigma$ STAN 2.84 <0.01 0.37 2.11 2.60 2.85 3.10 3.57
BPA 2.84 0.10 0.66 2.00 2.50 2.79 3.06 3.79

The computational time of the BPA to sample from the marginal posterior
distributions of both $\mu$ and $\log\sigma$ was measured to be 22.08 seconds,
which was 55% lower than the computational time of NUTS.

We note that there was a bias of 0.4 between the posterior mean estimates for
$\mu$ obtained from STAN and the BPA. Although this bias is small compared to
its magnitude (2% of the mean), this difference still needs to be accounted
for. This deviation can be explained by the arbitrary nature of the limits of
the uniform distribution as the prior. Simulating from a uniform distribution
with a different width affects the mean and the quantile estimates because of
the dynamics of the balls from points far away from the expected posterior
value. For this example, we used starting points simulated from $Unif(0,100)$,
but changing the upper limit to 50 or 1000 gave us different results. The
approximate posterior obtained from the BPA is slightly narrower compared to
the resulting posterior from NUTS. However, the estimated posterior median
from the BPA was consistent with the other estimation methods such as the
Random Walk Metropolis-Hastings MCMC, Independence sampling, and Gibbs
Sampling as calculated by Albert (2009).

The posterior distribution for the logarithm of the scale parameter, which we
refer to as the log-scale parameter, obtained from the BPA is similar to the
resulting posterior distribution from NUTS as shown in Figure 7. We also
observed a slightly heavier tail for the BPA, which after further
investigation was determined to be caused by the truncation set when we
simulated from the uniform prior. To account for this, we investigated the
posterior median instead of the posterior mean in comparing the two methods.
The posterior medians suggest that the BPA is shifted to the left compared to
the NUTS distribution. Based on the values provided by Albert (2009), the
posterior median resulting from NUTS was close to other Monte Carlo methods
(2.85) and other samplers such as the Gibbs sampler (2.86) while the posterior
median for the logarithm of the scale parameter obtained by the BPA (2.79) is
close to the same posterior median obtained when a Laplace normal
approximation (2.77) is used (Albert, 2009). However, the 95% credible
intervals are wider for the BPA approximation compared to the Laplace normal
approximation.

###  4.4 Discussion

The resulting distributions from the BPA are comparable to the distributions
obtained from the HMC implementation through the NUTS algorithm; however, a
notable difference between the results of the BPA and NUTS was that the BPA
generally yielded a narrower distribution as attested by the lower standard
deviations and quantile estimates. This trend was observed for distributions
with unbounded or half-bounded supports for its parameters such as the Poisson
distribution ($\lambda\in[0,\infty)$) and the Cauchy distribution
($\mu\in(-\infty,\infty)$ and $\sigma\in[0,\infty)$). It is imperative to know
where to truncate the improper priors when sampling for initial points using
the BPA since the resulting distributions are affected by where the chain
starting points are located. If the chosen starting point is far away from the
posterior mode, the chain might not converge. Aside from the nature of the
parameter supports, the choice of prior can also lead to narrower posterior
distributions after implementing the BPA. If the starting points are confined
within a narrow region of the parameter space, it might not be able to explore
regions that was covered by the HMC.

The BPA also reduced the computational time by at least 95% compared to the
rstan implementation for the single parameter models and 55% for the Cauchy
model. As stated in Section 3.2, the HMC entails simulating momentum values at
the start of each iteration to keep the value of the Hamiltonian constant.
With the BPA, there is no need to simulate velocity values after each
iteration since the trajectory is dictated by the probability of taking a
certain path given by the path integral formulation. Aside from the rejection
of the candidate point, the only instances where the velocities are resampled
are when the chain is stuck or the likelihood at the candidate point is
undefined. The fast convergence of the BPA also allows us to use shorter
chains which also leads to faster sampling from the target distribution. This
significant decrease in computational time shows the potential of the BPA for
sampling target distributions more efficiently compared to HMC and other
standard MCMC algorithms.

##  5 Conclusion and Next Steps

Using path integrals in the Bayesian analysis of Markov chains provided
valuable information on how the prior distribution evolves to the posterior
distribution. We found that the path followed by the Markov chain in the
parameter space is dictated by the value of the log-likelihood and the noise
strength of the Markov process. Using this information, we devised the BPA,
which is a Monte Carlo algorithm that utilizes a Lagrangian framework to
predict the trajectory of a Markov chain in the parameter space. We then
compared the BPA to the Hamiltonian Monte Carlo, which is another standard
sampling method that uses concepts from physical systems, as implemented by
the NUTS algorithm in the R package rstan. The BPA was found to be at least
95% faster than the NUTS algorithm in sampling posterior distributions for
single parameter models and 55% faster in sampling the marginal posterior
distributions for the Cauchy likelihood parameters. The resulting approximate
posterior distributions from the BPA is comparable to the ones generated using
NUTS, which is supported by visual and numerical evidences. Furthermore, the
standard deviation for the BPA posterior distributions are generally lower
than the NUTS posterior distributions.

The BPA was also applied to a multi-parameter Cauchy model to estimate the
location and the logarithm of the scale parameters. A slight bias was observed
for the estimated posterior means for the location parameter, but the
posterior median of the location parameter was found to be consistent between
other MCMC methods. As for the log-scale parameter, the estimated posterior
median from the result of the BPA is close to the posterior median when the
Laplace normal approximation is used. The distribution appears to be slightly
shifted to the left and is slightly wider than the posterior distribution
obtained from NUTS.

We aim to extend the BPA to sample from joint distributions for multi-
parameter models by accounting for multiple Markovian processes in the
parameter space. Sampling from the joint distribution of the parameters would
only need one implementation of the BPA instead of its current state where it
needs multiple implementations to sample from marginal distributions. We
believe that this will improve the performance of the algorithm in
approximating the exact posterior distributions of multi-parameter models.

## References

  * Albert (2009) Albert, J. (2009), Bayesian computation with R, Springer Science & Business Media.
  * Albert (2018) — (2018), LearnBayes: Functions for Learning Bayesian Inference, URL https://CRAN.R-project.org/package=LearnBayes. R package version 2.15.1.
  * Albeverio et al. (1976) Albeverio, S., Höegh-Krohn, R., and Mazzucchi, S. (1976), Mathematical theory of Feynman path integrals, vol. 523, Springer.
  * Andrieu et al. (2009) Andrieu, C., Doucet, A., and Holenstein, R. (2009), “Particle markov chain monte carlo for efficient numerical simulation,” in Monte Carlo and quasi-Monte Carlo methods 2008, Springer, pp. 45–60.
  * Andrieu et al. (2010) — (2010), “Particle markov chain monte carlo methods,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), 72, 269–342.
  * Betancourt (2017) Betancourt, M. (2017), “A conceptual introduction to hamiltonian monte carlo,” arXiv preprint arXiv:1701.02434.
  * Caldeira and Leggett (1983) Caldeira, A. O. and Leggett, A. J. (1983), “Path integral approach to quantum brownian motion,” Physica A: Statistical mechanics and its Applications, 121, 587–616.
  * Casella and George (1992) Casella, G. and George, E. I. (1992), “Explaining the gibbs sampler,” The American Statistician, 46, 167–174.
  * Chang et al. (2015) Chang, J. C., Fok, P.-W., and Chou, T. (2015), “Bayesian uncertainty quantification for bond energies and mobilities using path integral analysis,” Biophysical journal, 109, 966–974.
  * Chang et al. (2014) Chang, J. C., Savage, V. M., and Chou, T. (2014), “A path-integral approach to bayesian inference for inverse problems using the semiclassical approximation,” Journal of Statistical Physics, 157, 582–602.
  * Chen et al. (2003) Chen, Z. et al. (2003), “Bayesian filtering: From kalman filters to particle filters, and beyond,” Statistics, 182, 1–69.
  * Del Moral (1996) Del Moral, P. (1996), “Nonlinear filtering using random particles,” Theory of Probability & Its Applications, 40, 690–701.
  * Doucet et al. (2001) Doucet, A., De Freitas, N., and Gordon, N. (2001), “An introduction to sequential monte carlo methods,” in Sequential Monte Carlo methods in practice, Springer, pp. 3–14.
  * Duane et al. (1987) Duane, S., Kennedy, A. D., Pendleton, B. J., and Roweth, D. (1987), “Hybrid monte carlo,” Physics letters B, 195, 216–222.
  * Feynman et al. (2010) Feynman, R. P., Hibbs, A. R., and Styer, D. F. (2010), Quantum mechanics and path integrals, Courier Corporation.
  * Fujii and Hatakenaka (2018) Fujii, T. and Hatakenaka, N. (2018), “A path integral approach to bayesian inference in markov processes,” Global Journal of Pure and Applied Mathematics, 14, 721–731.
  * Gelman et al. (2013) Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2013), Bayesian data analysis, CRC press.
  * Geman and Geman (1984) Geman, S. and Geman, D. (1984), “Stochastic relaxation, gibbs distributions, and the bayesian restoration of images,” IEEE Transactions on pattern analysis and machine intelligence, 721–741.
  * Girolami and Calderhead (2011) Girolami, M. and Calderhead, B. (2011), “Riemann manifold langevin and hamiltonian monte carlo methods,” Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73, 123–214.
  * Glimm and Jaffe (2012) Glimm, J. and Jaffe, A. (2012), Quantum physics: a functional integral point of view, Springer Science & Business Media.
  * Goldstein et al. (2002) Goldstein, H., Poole, C., and Safko, J. (2002), “Classical mechanics,” .
  * Hawking (1979) Hawking, S. W. (1979), “The path-integral approach to quantum gravity,” in General relativity.
  * Hoffman et al. (2014) Hoffman, M. D., Gelman, A., et al. (2014), “The no-u-turn sampler: adaptively setting path lengths in hamiltonian monte carlo.” J. Mach. Learn. Res., 15, 1593–1623.
  * Holenstein (2009) Holenstein, R. (2009), “Particle markov chain monte carlo,” PhD thesis, University of British Columbia.
  * Lan et al. (2015) Lan, S., Stathopoulos, V., Shahbaba, B., and Girolami, M. (2015), “Markov chain monte carlo from lagrangian dynamics,” Journal of Computational and Graphical Statistics, 24, 357–378.
  * Liu and Chen (1998) Liu, J. S. and Chen, R. (1998), “Sequential monte carlo methods for dynamic systems,” Journal of the American statistical association, 93, 1032–1044.
  * MacKenzie (2000) MacKenzie, R. (2000), “Path integral methods and applications,” arXiv preprint quant-ph/0004090.
  * Pompe et al. (2020) Pompe, E., Holmes, C., and Łatuszyński, K. (2020), “A framework for adaptive mcmc targeting multimodal distributions,” The Annals of Statistics, 48, 2930–2952.
  * R Core Team (2017) R Core Team (2017), R: A Language and Environment for Statistical Computing, R Foundation for Statistical Computing, Vienna, Austria, URL https://www.R-project.org/.
  * Robert (2007) Robert, C. (2007), The Bayesian choice: from decision-theoretic foundations to computational implementation, Springer Science & Business Media.
  * Sakurai and Napolitano (2011) Sakurai, J. and Napolitano, J. (2011), “Modern quantum mechanics 2nd ed,” .
  * Smith and Roberts (1993) Smith, A. F. and Roberts, G. O. (1993), “Bayesian computation via the gibbs sampler and related markov chain monte carlo methods,” Journal of the Royal Statistical Society: Series B (Methodological), 55, 3–23.
  * Stan Development Team (2020) Stan Development Team (2020), “RStan: the R interface to Stan,” URL http://mc-stan.org/. R package version 2.21.2.
  * Szabados (2010) Szabados, T. (2010), “An elementary introduction to the Wiener process and stochastic integrals,” arXiv:1008.1510 [math], URL http://arxiv.org/abs/1008.1510. ArXiv: 1008.1510.
  * Tipler and Mosca (2007) Tipler, P. A. and Mosca, G. (2007), Physics for scientists and engineers, Macmillan.
  * Wiegel (1975) Wiegel, F. W. (1975), “Path integral methods in statistical mechanics,” Physics Reports, 16, 57–114.
  * Wio (2013) Wio, H. S. (2013), Path integrals for stochastic processes: An introduction, World Scientific.
  * Yuan et al. (2012) Yuan, K., Girolami, M., and Niranjan, M. (2012), “Markov chain monte carlo methods for state-space models with point process observations,” Neural Computation, 24, 1462–1486.