arXiv:1710.07755v1 [math.ST] 21 Oct 2017

A path integral approach to Bayesian inference
in Markov processes
Toshiyuki Fujiia , Noriyuki Hatakenakab,∗
a Department

of Physics, Asahikawa Medical University, Midorigaoka-higashi, Asahikawa,
078-8510, Japan
b Graduate School of Integrated Arts and Sciences, Hiroshima University, Higashi-Hiroshima
739-8521, Japan

Abstract
We formulate Bayesian updates in Markov processes by means of path integral
techniques and derive the imaginary-time Schrödinger equation with likelihood
to direct the inference incorporated as a potential for the posterior probability
distribution.
Keywords: Bayesian inference, Path integral technique, Imaginary-time
Schrödinger equation
2010 MSC: 62F15, 58D30, 60H15

1. Introduction
Statistical inference is a fundamental technology indispensable in modern
society. In particular, Bayesian inference is explosive growth in a wide range of
all kind of human activities, including philosophy, medicine, sport, law as well
as science and technology [1, 2]. This uses a procedure to update the probability
sequentially based on the Bayes’ theorem that incorporates new evidence into
prior probabilities to form updated probability estimates. The form of this
theorem is
P (H|E) =

P (E|H)
P (H),
P (E)

(1)

∗ Corresponding

author
Email address: noriyuki@hiroshima-u.ac.jp (Noriyuki Hatakenaka)

Preprint submitted to Elsevier

October 24, 2017

where P (E|H) is the probability for the evidence E given a hypothesis H,
P (H) is the prior probability that the hypothesis is true, P (H|E) is the poste-

5

rior probability that the hypothesis is true given the evidence, and P (E) =
P
H P (E|H)P (H) is the probability for the evidence averaged over all hy-

potheses. In this way, the Bayes’ theorem relates the prior probability of the

hypothesis before getting the evidence P (H) to the posterior probability of
the hypothesis after getting the evidence P (H|E) through the likelihood ratio
P (E|H)/P (E).
As is well known, the time evolution of the conditional probability (the
transition probability) of the Markov process follows the Chapman-Kolmogorov
equation [3, 4];
p(xi+1 , ti+1 |xi−1 , tt−1 ) =
10

Z

p(xi+1 , ti+1 |xi , ti )p(xi , tt |xi−1 , ti−1 )dxi ,

(2)

where p(xi , ti |xj , tj ) is the conditional probability density function defined as
the p(xi , ti |xj , tj )dxi dxj ≡ P (xi < Xi ≤ xi + dxi |xj < Xj ≤ xj + dxj ) with
xj and Xj being the realization and the random variable at a certain time
tj . This expresses the chain condition of the Markov process that the state
in the near future depends only on the state at the present time, not on past

15

memories. This is consistent with the idea of infinitesimal time evolution in path
integral techniques of quantum mechanics. In fact, Onsager and Machlup later
incorporated Markov chain conditions into path integral techniques [5]. Along
this line, this paper extends the Onsager-Machlup theory by incorporating the
likelihood required by Bayesian inference into the Markov chain conditions.

20

2. A path integral representation of Bayesian updates
Suppose the stochastic processes under the probability distribution density
functions pi (qi ) and pi (xi ) that are defined as pi (qi )dqi ≡ P (qi < Qi ≤ qi + dqi )
and pi (xi )dxi ≡ P (xi < Xi ≤ xi + dxi ), respectively. Here the continuous
random variables Qi and Xi act as the roles of H and E in the Bayes’ theorem

2

in the i-th update. Starting with the law of total probability
Z
pi+1 (qi+1 ) = dqi p(qi+1 |qi )pi (qi ),

(3)

the first Bayesian update is expressed as
Z
p(x1 |q1 )
p1 (q1 |x1 ) = dq0
p(q1 |q0 )p0 (q0 ).
p(x1 )

(4)

After repeating Bayesian updates of N times, the posterior probability distribution density function is given as
Z
pN (qN |x1 , x2 , · · · xN ) = dq0 K(qN , q0 ; x1 , x2 , · · · xN )p0 (q0 ),

(5)

where the kernel K(qN , q0 ; x1 , x2 , · · · xN ) is given as
K(qN , q0 ; x1 , x2 , · · · xN )
Z
Z
Z
N
Y
p(xi |qi )
= dq1 dq2 · · · dqN −1
p(qi |qi−1 ).
p(xi )
i=1

(6)

This describes repeated Bayesian updates representing the time evolution of
stochastic processes starting from the prior probability distribution density function p0 (q0 ) inferred from the repeated conditional probabilities. Hereafter, we
simply represent pi (qi |x1 , x2 , · · · xi−1 ) = pi (qi |x) where x denotes the set of the
25

realizations.
Here, we introduce the parameter ti = iǫ with a small quantity ǫ and an
integer i, in order to prepare for converting discrete time to continuous time
such as qi = q(ti ). In this notation, for example, the probability distribution
density function and the transition probability distribution density function are

30

expressed as pi (qi = q|x) = p(q, ti |x) and p(qj = q ′ |qi = q) = p(q ′ , tj |q, ti ),

respectively. Let us consider the infinitesimal time evolution of p(q ′ , tj |q, ti )

from time ti to tj = ti+1 . Since p(q ′ , ti+1 |q, ti ) = p(q ′ , ti + ǫ|q, ti ), then
p(q ′ , ti + ǫ|q, ti )

= p(q ′ , ti |q, ti ) + ǫ

∂p(q ′ , t′i |q, ti )
∂t′i

+ O(ǫ2 ).

(7)

t′i =ti

Suppose that the stochastic process that we are assuming is a Wiener process.
If the stochastic process is a Markov process, the random variable q obeys the
3

Ito-type stochastic partial differential equation and the probability distribution
density function based on this random variable follows the Fokker-Planck (FP)
equation. Furthermore, the FP equation results in a diffusion equation in the
case of a Wiener process leading to Brownian motion which has no driving force.
Therefore, the second term on the right hand side of Eq. (7) is replaced by the
diffusion equation;
∂p(q ′ , t′ |q, t)
1 ∂ 2 p(q ′ , t′ |q, t)
=
,
D
∂t′
2
∂q ′2

(8)

where D is the diffusion constant. The resulting infinitesimal time-evolution
equation is
p(q ′ , ti + ǫ|q, ti ) =



ǫ ∂2
1 + D ′2 p(q ′ , t′i |q, ti )
2 ∂q

+ O(ǫ2 ).

(9)

t′i =ti

By using the relation p(q ′ , t′i |q, ti )|t′i =ti = δ(q ′ − q), this is reduced to


ǫ ∂2
′
p(q , ti + ǫ|q, ti ) = 1 + D ′2 δ(q ′ − q) + O(ǫ2 )
2 ∂q
 Z


ǫ ∂2
ds i(q′ −q)s
+ O(ǫ2 )
= 1 + D ′2
e
2 ∂q
2π
Z

′
ds 
ǫ
=
1 + D(−s2 ) ei(q −q)s + O(ǫ2 )
2π
2
Z
ds − ǫ Ds2 i(q′ −q)s
e
+ O(ǫ2 )
e 2
=
2π
(q′ −q)2
1
=√
e− 2Dǫ + O(ǫ2 ).
2πDǫ

(10)

Here we define the velocity of the random variable by q̇ = (q ′ − q)/ǫ. The
equation is then given as
q̇2
1
p(q ′ , ti + ǫ|q, ti ) = √
e− 2D ǫ + O(ǫ2 ).
2πDǫ

(11)

This is consistent with the Onsager-Machlup formula [5] in the path integral
representation of stochastic processes.
Now, let us incorporate the likelihood which is the most significant factor in
Bayesian update into the path integral form. We adopt the exponential function
used in the previous studies [6] as likelihood p(x, ti |q, ti );
p(x, ti |q, ti ) = M exp[−V (q(ti ) − x(ti ))ǫ],
4

(12)

where M is a normalization constant. This form satisfies the requirement for
the probability density function that it has a non-negative value in the range
from −∞ to ∞. Thus, the kernel is represented as
K(q(tN ), q(t0 ), {x(ti ), i = 1, 2, · · · N })
Z
Z
Z
N
=M
dq(t1 ) dq(t2 ) · · · dq(tN −1 )
"
#
N 
X
q̇(tn )2
× exp −ǫ
+ V (q(tn ) − x(tn ))
.
2D
i=1

(13)

Now let us take a limit ǫ → 0 under the condition tN = N ǫ ≡ t′ . The discrete
variable ti = iǫ becomes dense at time interval [0, t′ ]. The set of realizations in

stochastic processes q and x can be represented by continuous functions q(t) and
x(t). Thus, the probability distribution density function pi (qi ) is also expressed
as p(q, t). In this limit. Eq. (13) is just a path integral description of Bayesian
updates;
K(q(t′ ), q(t), x(t)) =
with

D[q]e−S[q(t),q̇(t),x(t)] ,

(14)


q̇(t)2
+ V (q(t) − x(t)) dt,
(15)
S[q(t), q̇(t), x(t)] =
2D
0
R
R
R
R
where D[q] ≡ M N dq(t1 ) dq(t2 ) · · · dq(tN −1 ). This action S describes
Z

35

t′

Z



a particle moving in the random potential V . This is a main result of this

paper. The likelihood, which is the most important factor in Bayesian update,
is incorporated into the path integral form as a potential. This result seems to be
quite natural. As is well known, the potential leads the force, which determines
40

the direction of particle movement in mechanics. It is exactly consistent with
the role of likelihood in Bayesian update.

3. An imaginary-time Schödinger equation
Let us derive the dynamical equation of the posterior probability distribution
density associated with the Bayesian update based on the path integral formula
45

for Bayesian updates derived in the previous section. Under an infinitesimal

5

interval ǫ from t to t′ ( t′ = t + ǫ),
′

′

p(q , t |x) = M

Z

dqe

−ǫ



1
2D

n

q′ −q
ǫ

o2

+V (q−x(t))



p(q, t|x).

(16)

Let us take the Taylor series expansions for probability distribution density
functions on both sides of the equation. The probability distribution density
function p(q ′ , t′ |x) is expanded with respect to ǫ as
p(q ′ , t + ǫ|x) = p(q ′ , t|x) +

∂p(q ′ , t|x)
ǫ + O(ǫ2 ).
∂t

On the other hand, the probability distribution density function p(q, t|x) on the
right hand side can be expressed as p(q ′ − η, t|x) to meet the variable on the left
side by introducing a variable η = q ′ − q, and is expanded with respect to η,
p(q ′ − η, t|x)
= p(q ′ , t|x) −

1 ∂ 2 p(q ′ , t|x) 2
∂p(q ′ , t|x)
η+
η + O(η 3 ).
′
∂q
2
∂q ′2

Then, we obtain
∂p(q ′ , t|x)
p(q ′ , t|x) +
ǫ
∂t
Z
η2
= M dηe− 2ǫ [1 − ǫV (q − x(t))]


∂p(q ′ , t|x)
1 ∂ 2 p(q ′ , t|x) 2
′
× p(q , t|x) +
η+
η .
∂q ′
2
∂q ′2

(17)

After Gaussian integration with
Z ∞
√
η 2n exp(−η 2 /2Dǫ)dη = (2n − 1)!!(Dǫ)n 2πDǫ,
−∞

we get
p(q ′ , t|x) +

∂p(q ′ , t|x)
ǫ
∂t



√
D ∂ 2 P (q ′ , t|x)
′
ǫ
= M 2πDǫ [1 − ǫV (q − x(t))] p(q , t|x) +
2
∂q ′2
√
= M 2πDǫp(q ′ , t|x)


√
D ∂ 2 p(q ′ , t|x)
+ O(ǫ2 ).
+ M 2πDǫ3 −V (q − x(t)) +
2
∂q ′2
6

(18)

Comparing the first term on both sides of equation yields M = (2πDǫ)−1/2 .
The remaining terms eventually form


∂p(q, t|x)
D ∂2
−
+ V (q − x(t)) p(q, t|x).
= −
∂t
2 ∂q 2

(19)

This describes the time evolution of the posterior probability density function
in Bayesian inference in Markov processes. This is nothing but the Schödinger
equation in quantum mechanics except that it is imaginary time. This equa50

tion is very advantageous for Bayesian inference in a way different from the
path integral form because there are many technological accumulations in the
Schrödinger equation along with the development of quantum mechanics. In
the next section, we will look at some of these advantages.
3.1. Example: Gaussian likelihood
Let us consider the time evolution of the probability distribution density
function using the imaginary-time Schrödinger equation for Bayesian update
obtained in the previous section by adopting the Gaussian likelihood as a concrete example given by
V =

[q − x(t)]2
.
2∆2

(20)

The imaginary-time Schrödinger equation for Bayesian update to be solved is
then
−

∂p(q, t|x)
=
∂t



[q − x(t)]2
D ∂2
p(q, t|x).
+
−
2 ∂q 2
2∆2

(21)

This equation is the imaginary time version of the Schrödinger equations of
forced harmonic oscillators with a random force x(t)
F =−
55

dV
q
q
= − 2 q + 2 x(t).
dq
∆
∆

(22)

This partial differential equation can simply be solved using the Feynman-Kac
formula [7, 8]. However, we will seek the solution using the operator method
developed in quantum mechanics to show the usefulness of the imaginary-time
Schrödinger equation in Bayesian inference.

7

By introducing new operators called creation (a) and annihilation (a† ) operators in quantum mechanics;
s √
∆ D
q=
(a + a† ),
2
s
∂
1
√ (a − a† ),
=
∂q
2∆ D

(23)
(24)

with [a, a† ] ≡ aa† − a† a = 1, the equation is expressed by


s√
√

 x(t)2
∂p(q, t|x)  D †
D
†
−
a a + 1/2 +
+
x(t)
(a
+
a
)
=
p(q, t|x).

 ∆
∂t
2∆2
2∆3

(25)

The first two terms on the right hand side represent operators corresponding
to Hamiltonians for a simple harmonic oscillator and a random force in quantum mechanics. On the other hand, the third term expresses their interaction.
We focus on the time evolution of the probability distribution density function
p′ (q, t|x) due to the interaction using the operator
)
( √ 

Z t
1
x(t′ )2 ′
− D
†
a a+
t−
dt ,
U0 (t) = exp
2
∆
2
0 2∆
associated with p(q, t|x) = U0 p′ (q, t|x), this reads
s√
′
D
∂p (q, t|x)
−
U0 (−t)(a + a† )U0 (t)p′ (q, t|x).
= x(t)
∂t
2∆3

(26)

(27)

According to the Baker-Campbell-Hausdorff formula [9] for any two operators,
say A and B,
1
1
eA Be−A = B + [A, B] + [A, [A, B]] + · · · [A, [A, · · · [A, B] · · · ] + · · · , (28)
2
n!
the terms U0 (−t)(â + â† )U0 (t) are expressed as
U0 (−t)(a + a† )U0 (t) = e−

√

Dt
∆

a+e

√
Dt
∆

a† .

Thus, Eq. (27) is now represented as
s√

√
Dt
D  − √Dt
∂p′ (q, t|x)
∆ a + e ∆ a†
p′ (q, t|x)
e
= x(t)
−
∂t
2∆3
≡ H(t)p′ (q, t|x).
8

(29)

(30)

The formal solution of this equation is well known and given as


 Z t
′
′
′
p′ (q, 0),
p (q, t|x) = T exp −
H(t )dt

(31)

0

where T is the time ordering operator


A(t)B(t′ ) for t > t′
′
T [A(t), B(t )] =

B(t′ )A(t) for t < t′ .

This is also expressed by [10]
" Z
t

′

p (q, t|x) = exp −

1
H(t )dt −
2
′

0

′

Z

t

dt

0

′

Z

0

(32)

t′
′′

′

′′

#

dt [H(t ), H(t )] p(q, 0), (33)

since [H(t), [H(t′ ), H(t′′ )]] = 0 from
√ h
√
√
√
√
i
Dt′
Dt′
Dt′′
Dt′′
D
′
′′
′
′′
[H(t ), H(t )] = x(t )x(t ) 3 (e− ∆ a + e ∆ a† ), (e− ∆ a + e ∆ a† )
2∆
√ h
√
√
i
D(t′ −t′′ )
D(t′ −t′′ )
D
′
′′
∆
∆
= x(t )x(t ) 3 (e−
[a, a† ] + e
[a† , a]
2∆
!
√
√
D
D(t′ − t′′ )
′
′′
.
(34)
= −x(t )x(t ) 3 sinh
∆
∆
The time evolution of the probability distribution density function which focused
on the interaction is then given as

s√

Z t

√
√

′
′
Dt
Dt
D
p′ (q, t|x) = ef (t) exp −
dt′ x(t′ )
ae− ∆ + a† e ∆ p′ (q, 0), (35)
3
2∆
0
√
√
Rt
R t′
where f (t) ≡ ( D/2∆3 ) 0 dt′ 0 dt′′ x(t′ )x(t′′ ) sinh( D(t′ − t′′ )/∆). This can

be reverted to the posterior probability distribution density function using the
relation p(q, t|x) = U0 (t)p′ (q, t|x), resulting in
#
" Z
r
t
√
√
Dt′
Dt′
1
−1
p(q, t|x) = ef (t) exp −
dt′ x(t′ )
U0 (t)(ae− ∆ + a† e ∆ )U0
2∆3
0
× U0 (t)p(q, 0)


s√
Z t
√
√
′ −t)
′ −t)
D(t
D(t
D
+ a† e ∆ )
= ef (t) exp −
(ae− ∆
dt′ x(t′ )
3
2∆
0
× U0 (t)p(q, 0),

9

(36)

where the relations p(q, 0) = p′ (q, 0) and U0 (t)F (A)U0−1 (t) = F (U0 (t)AU0−1 )
60

that holds for any function F (A) of the operator A are used.
The prior probability distribution density function density is assumed to be
given by the superposition of the harmonic oscillator solutions fn as follows
p(q, 0) =

X

cn fn (q),

(37)

n

where
2
1
− q√
e 2∆ D Hn
fn (q) =  √ 1/4 √
2n n!
π∆ D

q
p √
∆ D

!

,

(38)

with Hn (x) being the n-th order Hermite polynomial at x. The expansion coefR
ficient cn is given by cn = dqfn (q)P0 (q, 0). Starting with this prior probabil-

ity distribution density function, the posterior probability distribution density
function is given as follows;


s√
Z t
√
√


D(t′ −t)
D(t′ −t)
D

+ a† e ∆
ae− ∆
p(q, t|x) = ef (t) exp −
dt′ x(t′ )
2∆3
0
!
∞
X
× U0 (t)
cn fn (q)
= ef (t)−

x(t′ )2
0 2∆2

Rt

n=0
∞
X
dt′

cn e

−

√
D
∆

(n+ 21 )t

n=0



× exp −

Z

0

t


s√
√
√


D(t′ −t)
D(t′ −t)
D
 fn (q).
+ a† e ∆
ae− ∆
dt′ x(t′ )
2∆3

(39)

Note that
e

−

Rt
0

dt′ x(t′ )

q√

D
2∆3

√
√


D(t′ −t)
D(t′ −t)
∆
∆
ae−
+a† e

10

= e−x

′′

∂
g(t)
(t)q −x′ (t) ∂q

e

e

,

(40)

where
√ Z t
o
n√
D
x (t) =
D(t − t′ )/∆ ,
dt′ x(t′ ) sinh
∆ 0
)
(√
Z t
D(t − t′ )
1
′′
′
′
x (t) = 2
,
dt x(t ) cosh
∆ 0
∆
√ Z t
Z t
D
′
g(t) =
dt
dt′′ x(t′ )x(t′′ )
2∆3 0
0
)
(√
)
(√
D(t − t′ )
D(t − t′′ )
sinh
,
× cosh
∆
∆
′

(41)
(42)

(43)

∂

and ec ∂q f (q) = f (q + c).
As a result, we reach the final expression for the posterior probability distribution density function p(q, t|x) as
p(q, t|x) = ef (t)+g(t)−
×

∞
X

Rt
0

x(t′ )2
2∆2

cn e −

n=0

dt′ −x′′ (t)q

√
D
∆

e

(n+ 21 )t f (q − x′ (t)).
n

(44)

We can find the role of Bayesian update in stochastic process in this solution
obtained from stationary solution fn of the Schrödinger equation. Evidence
incorporated in Bayesian update is embedded in the argument of the function
65

fn . This means that the center of the posterior distribution density function
is shifted with every update. In addition, the first n-independent exponent
is interpreted as the global phase corresponding to the geometric phase. These
features are difficult to extract in the path integral form like Feynman-Kac form.
Therefore, this is the greatest advantage of analysis that makes it possible to

70

elucidate elementary processes of stochastic processes using the Schrödinger
equation.

4. Summary
We have formulated the Bayesian inference in Markov processes by means
of path integral techniques and have succeeded in incorporating the likelihood
75

which is the key in Bayesian update as a potential into the path integral form.
11

Then we have derived the imaginary-time Schrödinger equation from path integral representation of the evolution of the Bayesian inference in Markov processes. Through an example of Gaussian likelihood, we have showed that the
operator technique for the Schrödinger equation is effective for elucidating the
80

elementary processes of stochastic processes.

Bibliography
[1] A. O’Hagan, M. West, The Oxford handbook of applied Bayesian analysis,
Oxford University Press, 2010.
[2] J. K. Kruschke, T. M. Liddell, The bayesian new statistics:
85

Hy-

pothesis testing, estimation, meta-analysis, and power analysis from
a bayesian perspective, Psychonomic Bulletin & Review (2017) 1–
29doi:10.3758/s13423-016-1221-4.
[3] S. Chapman, On the brownian displacements and thermal diffusion of
grains suspended in a non-uniform fluid, Proceedings of the Royal Society

90

of London A: Mathematical, Physical and Engineering Sciences 119 (781)
(1928) 34–54. doi:10.1098/rspa.1928.0082.
[4] A. Kolmogoroff, Über die analytischen methoden in der wahrscheinlichkeitsrechnung, Mathematische Annalen 104 (1) (1931) 415–458.
doi:10.1007/BF01457949.

95

[5] L. Onsager, S. Machlup, Fluctuations and irreversible processes, Phys. Rev.
91 (1953) 1505–1512. doi:10.1103/PhysRev.91.1505.
[6] W.

Bialek,

learning

C.

G.

probability

Callan,

S.

distributions,

P.

Strong,

Field

Phys.Rev.Lett

theories
(1996)

for

4693–

4697doi:10.1103/PhysRevLett.77.4693.
100

[7] M. Kac, On distributions of certain wiener functionals, Transactions of the American Mathematical Society 65 (1) (1949) 1–13.
doi:10.1090/S0002-9947-1949-0027960-X.
12

[8] B. K. Cheng, The exact propagator of time-dependent forced harmonic
oscillator, Revista Brasileira cie Ffsica 13 (1) (1983) 220–230.
105

[9] H. Baker, Alternants and continuous groups, Proceedings of the London
Mathematical Society 2 (1) (1905) 24–47. doi:10.1112/plms/s2-3.1.24.
[10] E. Merzbacher, Quantum Mechanics, Wiley, 1997.

13

