Essential Questions for the Exam 2018, AMCS 308,

Stochastic Methods in Engineering

May 21, 2018

1. Formulate the basic properties a Wiener process.

The one dimensional Wiener process is a mapping W : [0, ∞) × Ω → R, where Ω
is a probability space with probability measure associated. This Process satisfies the
following properties
• W0 = 0.
• The mapping t → Wt is almost surely continuos in [0, ∞).
• The increments of Wt are independent and normal distributed random variables:
Wt − Ws ∼ N (0, t − s)f or0 ≤ s < t.
2. State and derive the Central Limit Theorem. State the BerryEsseen The-
orem.
The Central
PLimit Theorem: Let X1 , . . . , Xn be IID with mean µ and variance σ 2 . Let
−1 n
Xn = n i=1 Xi . Then
√
Xn − µ n(X n − µ)
Zn = q = →Z
σ
V [X n ]

where Z ∼ N (0, 1).

Proof : Let Yi = (Xiσ−µ) . Then, Zn = n−1/2 ni=1 Yi . Now
P
n
Pn if φYi (t) is the moment-
generating
√ function of Yi , then (φYi (t)) is the MGF of i=1 Yi and then, φZn (t) =
(φYi (t/ n))n is the MGF of Zn . Now I have that:
(
φ0Yi (0) = E[Yi ] = 0,
φ00Yi (0) = E[Yi2 ] = V[Yi ] = 1.

Then:
t2 00 t3
φYi (t) = φYi (0) + tφ0Yi (0) + φYi (0) + φ000 (0) + . . .
2! 3! Yi
t2 t3 000
= 1 + 0 + + φYi (0) + . . .
2! 3!
t2 t3 000
= 1 + + φYi (0) + . . . ,
2! 3!
1
 and then:
 t n
φZn (t) = φYi √
n
2
 t t3 000
n 2
= 1+ + 3/2
φYi (0) + . . . → et /2 ,
2!n 3!n
which is the MGF of a N (0, 1).
The Berry-Essen Inequality:
let’s define,

E[|Y − EY |3 ]
λ= <∞ (1)
σy
then the inequality is given by

c∗ λ3
|F (x) − φ(x)| ≤ √ (2)
(1 + |x|)3 M

where φ is the CDF of a standard normal and c∗ is a constant.

3. Show how to apply the Monte Carlo method to compute an integral and
discuss the corresponding error.
Consider,
Z 1 M
1 X
f (x) dx ≈ f (xi ) (3)
0 M i=1
where xi ∼ U [0, 1] IID samples.
The variance is given by,
M M
1 X 1 X V [f (xi )]
V[ f (xi )] = 2 Vf (xi ) = (4)
M i=1 M i=1 M

4. Describe and give all the proofs on the sampling method of acceptance-
rejection.
The method is as follows. We choose a proposal density η which has the share the
support of our target distribution f and covers it. Then we sample U ∼ U [0, 1]. And
for any  ∈ [0, 1], ( f (x)
 η(x) < U Reject the sample
f (x) (5)
 η(x) ≥ U Accept the sample
Proof:
f (x)
Z Z  η(x) Z Z
f (xk ) f (x)
P (Uk ≤  )= dµ dx =  η(x) dx =  f (x) dx =  (6)
η(xk ) 0 η(x)

2
 Let xk be the first accepted sample. We would like to show that xk is from the density
f . Consider,
k−1
X f (xk ) Y f (xm )
P (xk ∈ B) = P (xk ∈ B, Uk ≤  ) P (Um >  )
k≥1
η(xk ) m=1 η(xm )
| {z }
=1−
(7)
f (xk ) X
= P (xk ∈ B, Uk ≤  ) (1 − )k−1 ∵ IID
η(xk ) k≥1
| {z }
1/

And we have,
Z Z f (x) Z
f (xk ) η(x) f (x)
P (xk ∈ B, Uk ≤  )= dµ η(x) dx =  η(x) (8)
η(xk ) B 0 B η(x)
Hence, Z
f (xk )
P (xk ∈ B|Uk ≤  )= f (x) dx (9)
η(xk ) B

And so the sample are from the density f given that we have accepted them.

We look at the expected number of sampling from the proposal before we accept one
sample,
X X 1
E[K] = kP (K = k) = k(1 − )k−1  = (10)
k≥1 k≥1


We see that, ( →0

E[K] −−→ ∞
→1 (11)
E[K] −−→ 1

5. Motivate the use of variance reduction techniques. Describe two techniques

and discuss their computational efficiency.
Monte Carlo’s error is given by ,
M
r
 1 X  VY
E[Y ] − Yj  ≈ C α (12)

M j=1 M

We would like to reduce the variance without introducing bias.

(a) Control Variates: We would like to estimate EY . We sample Y and another

auxiliary RV, X for which EX is given and X and y are highly correlated.
For a given β > 0, we consider the unbiased estimator
M M
1 X β X
Vm = Yj − (Xj − EX) (13)
M j=1 M j=1

3
 Table 1: Comparison between MC and MC with Control Variates
Method #Samples Computational Work
c 2 c
MC ( T OL ) VY ( T OL )2 VY
Control Variates ( T OL )2 V[Vm (β ∗ )]
c c
(1 + ρ)( T OL )2 VY

We observe that,
1 1 
V[Vm ](β) = V[Y − βX] = VY + β 2 VX − 2βCov[X, Y ] (14)
M M
We find the minimizing β ,
Cov(Y, X)
β∗ = (15)
VX
Then,
∗ 1  Cov(Y, X)  1
V[Vm ](β ) = V[Y ] 1 − ≤ V[Y ] (16)
M V [X]V[Y ] M
Assume that the work to sample (X, Y ) is (1 + ρ) times the work of sampling Y
alone.
(b) Antithetic variates: Let Y = g(x) such that X has a symmetric distribution and
so EX = 0. Then we have X and −X identically distributed, which means that
Eg(X) = Eg(−X) and
h g(X) + g(−X) i
EY = E (17)
2
then we have the estimator,
M
1 X g(Xj ) + g(−Xj )
Vm = (18)
M j=1 2

Which has variance,

V[g(x) + g(−X)]
V[Vm ] = (19)
4M
Assume that the extra evaluation of g(−X) has the same cost as evaluating g(X).
Then we require that,

2V[Vm ] ≤ V[X] =⇒ V[g(X) + g(−X)] ≤ 2V[g(X)] (20)

and it sufficient that,

Cov[g(X), g(−X)] < 0 (21)
For linear g we have V[g(X) + g(−X)] = 0. As,

V [g(X) + g(−X)] = V [g(X) − g(X)] = 2V[g(x)] + 2 Cov(g(X), −g(X)) = 0 (22)

| {z }
=−V[g(X)]

4
6. Define the Poisson counting process and explain how to sample it.
It is a collection of random variables where N (t) is a counting of events that has
occurred upto time t starting from t0 . The number of events between time t and time
s is given by N (s) − N (t) which follows a poisson distribution.
Def.
(a) N (0) = 0
(b) N has independent increments
(c) If I = (s, t) ⊆ R for k = 1, 2, . . .
(λ|I|)k
P (N (t) − N (s) = k) = P (k jumps in I) = e−λ|I| (23)
k!
Where λ > 0 is the intensity of the process.
Given a sequence time t where t0 = 0 and tn = T with λ > 0:
(a) Sample ∆Ntn ∼ P oisson(λ(tn − tn−1 )) for k = 1, 2, . . . N
P
(b) Nt = ∆Ntk
k=1

7. Define the compound Poisson process and explain how to sample it.

A compound poisson process with rate λ > 0 and jump size distribution f is a contin-
uous time stochastic process {X(t)} given by
N (t)
X
X(t) = Zi
i=1

Where N (t) is a poisson counting process with rate λ > 0. The random variables Z
are IID with distribution f .
To sample a compound poisson process with jump amplitudes N (0, 1) and with parti-
tions in time ∆tn :
(a) Sample ∆Ntn ∼ P oisson(λ(tn − tn−1 )) ∀n = 1, 2, . . . N .
(b) Sample ∆Xtn ∼ N (0, ∆Nt2n ) ∀n = 1, 2, . . . N
N
P
(c) X(t) = ∆Xtn
n=1

8. What is a Brownian bridge? Discuss an algorithm to sample the mid point

of the Brownian bridge.

Given a Wiener process

 W (t) on [t1 , t2 ]. Let the Brownian bridge have the conditional
value B(s) = W (s) for s ∈ [t1 , t2 ]. Given B(t1 ) = a and B(t2 ) = b, the

W (t2 )=b
distribution of B at t ∈ [t1 , t2 ] is Normal with parameters,
t − t1 (t2 − t)(t − t1 )
E[B(t)] = a + (b − a) V[B(t)] = (24)
t2 − t1 (t2 − t1 )

5
9. Describe Kernel density estimates for non parametric estimation. Discuss
the choice of the window parameter h and the resulting error in the esti-
mate of the pdf. How does this discussion depend on the dimension of the
random variable?

We want to estimate the density ρy (y) given a set {Yl }M

l=1 IDD observations from a
common density ρy . The kernel density estimator of ρy of the point y is given by,
M y − Y 
1 X j
ρ̂k (y) = Ker (25)
M k j=1 k

if Ker(·) is a pdf, i.e. ( Ker(z) ≥ 0 and normalized). Then ρ̂k (y) will also be a candi-
date pdf estimator.
Next, we will show that no unbiased estimator can exit for all continuous ρy (y),

We look for a candidate with the least mean square error. An asymptotically unbiased
candidate estimator ρ̂km (y) such that,
M →∞
E[ρ̂km (y)] −−−−→ ρ(y) ∀y ∈ R (26)

We choose the window size parameter k > 0 to minimize the total error in the estimate
of Y ∼ ρ, given by

ρ̂k (y) − ρ(y) = E[ρ̂k (y)] − ρ(y) + ρ̂k (y) − E[ρ̂k (y)] (27)
| {z } | {z }
Bias sampling error

We look closer at the bias error,

Z
E[ρ̂k (y)] − ρ(y) = Kerk (y − z)ρ(z) dz − ρ(y)

= (Kerk ∗ ρ)(y) − ρ(y)

Z (28)
= Kerk (z)ρ(y − z) dz − ρ(y)
1 00
≈ σk2 k 2 ρ (y)
2
Note that we have expanded ρ(y − z) using a Taylor series. Where the first term can-
cels out with ρ(y) as the kernel is a pdf and the second term
R 2 is zero as the kernel is
2
symmetric.. And so we write the third term with σKer = z Kerk (z) dz plus higher
order terms.

Next we look at the sampling error,

By chebychev’s inequality, for any a > 1,
  1
P |EX − X| ≥ aσ(X) ≤ 2 (29)
a
6
We can approximate,  
E[ρˆk ](y) − ρ̂k (y) ≈ Cσρ̂k (y) (30)
 

And,
1 y − z 
V[ρ̂k ](y) = V[Ker ]
M k2 k
Z 0 Z 00 Z
ρ(y) 2 ρ (y) 2 ρ (y)
= Kerk (z) dz + zKer (z) dz + z 2 Ker2 (z) dz + . . .
Mk M 2M
(31)

Where we have expanded ρ(y − kz) using Taylor around kz.

Next we choose k that minimizes the total squared error,
Z
 2 1 4 00 2 ρ(y)
E[ρ̂k ](y) − ρ(y) + CV[ρ̂k ](y) = σKer (ρ (y)) + C Kerk2 (z) dz + . . . (32)
4 MK
Hence, k ∗ ∝ M −1/5 . And so the approximation error ∝ M −2/5 .However, we see that
this choice depends on y. Next we find an optimal k independent of y. This can be
done by minimizing the L2 error,
Z Z Z
2 2
E[(ρ(y) − ρ̂k (y)) ] dy ≤ 2 (ρ(y) − E[ρ̂k ](y)) dy + 2 V[ρ̂k ](y)
Z Z (33)
1 4 4 00 2 2c 2
using chebychev inq. ≤ σKer k |ρ (y)| dy + Kerk (z) dz + . . .
2 MK

We consider the Mean Square Error,

Z
M SE(ρ̂k ) = E (ρ̂k (y) − ρ(y))2 dy
Z
= E(ρ̂k (y) − ρ(y))2 dy
Z
= M SE(ρ̂k (y)) dy
Z Z Z Z
1 1
= Ker (y) dy + (1 − ) (ρ̂k ∗ ρ) (y) − 2 (Kerk ∗ ρ)(y)ρ(y) dy + ρ2 (y) dy
2 2
Mk M
1 1 00 1
= R(Ker) + k 2 σKer
2
R(ρ ) + O( + k4)
Mk 4 Mk
(34)
R 00 R 00
where R(Ker) = Ker2 (y) dy and R(ρ ) = (ρ )2 (y) dy. Hence,
1 1 00
M SE(ρ̂k ) ≈ R(Ker) + k 2 σKer
2
R(ρ ) (35)
Mk 4
Then the optimal choice of k is given by ,
R(Ker) 1/5 −1/5
kM SE = [ ] M (36)
4
σKer (ρ00 )

7
 The optimal value of k ∗ for a kde of an Rd random variable X, based on observing
X1 , X2 , . . . , Xn is as follows:

Let φ be a density. Then we have the estimate,

1 x
φH (x) = d
φ( ) (37)
H H
Set ,
n n
1X 1 X
fn = fX1 ,X2 ,...,Xk = φH (x − Xi ) = φ(H(x − Xi )) (38)
k i=1 k|H| i=1
where H is a smooth matrix to be specified. Also note that f is a pdf. Hence, the
MSE is given by ,

M SE(fˆn ) = Bias(fˆn )2 + Vfˆn

 1 Z  1
∗= φ(H −1 (y − x)f (x) dx − f (y)) + Vf [φ(H −1 (y − Xi )]
|H| k|H|2
1 2 1  1 Z 
≈ t t 2
Eφ [x H σ f (y)Hx] + 2 ˆ
φ (x)f (y − Hx) dx E[fn (y)]2
2 k |H|
(39)
Where we have used again Taylor expansion and a change of variables. We can proceed
to write,
1  2 f (y) Z ∇f (y)H 1 ˆ
M SE(fˆn ) = t 2 t
tr(H σ f (y)H)Eφ [xx ] + φ2 (x) dx − E[fn (y)]2
4 k|H| k|H| k
4
k f (y)
≈ [tr(∇2 f (y)Eφ [xxt ])]2 + d
||φ||22
4 kk
(40)
1 4
Hence, the optimal k ∗ ∝ M − d+k . And so the optimal choice ∝ M − d+k
10. What is affine prediction? Describe its use for the approximation of the
conditional expectation E[Y |Z]? When is this approximation exact?

The affine prediction is a method used to solve prediction problems, that means given
the joint distribution of a random variable (Y, Z) and the observations of Z, we want
to compute V̂ , the best prediction of Y .
In particular, consider the space:

Lp = {g : g is a polynomial with degre at most p} (41)

The affine prediction returns an Lp approximation of E[Y |Z]

.
Affine prediction: If we restrict ourselves to the L1 approximation of E[Y |Z]. So
lets consider, Y = g(Z) = at (Z − EZ) + b. To find (a∗ , b∗ ), we solve for m = 1, . . . , M .

8
 E[(at (Z − EZ) + b)(Zm − EZm )] = E[Y (Zm − EZm )] (42)
and
E[at (Z − EZ) + b] = E[Y ] (43)

This can be written in a matrix form as,

    
CovZ 0 a E[Y (Z − EZ)]
=
0 1 b EY

So we have that b∗ = EY and a∗ = Cov −1 (Z)E[Y (Z − EZ)] then,

g1∗ (Z) = E[Y (Z − EZ)]t Cov −1 (Z)(Z − EZ) + EY

 t (44)
= E[Y Z] − EY EZ Cov(Z)−1 (Z − EZ) + EY

11. Motivate the use of non-linear weighted least squares starting from a max-
imum posterior, Gaussian likelihood, Bayesian formulation.
The least squares technique for estimating M parameters X using N data D is given
by,
N
X
min X 2 = R2 (45)
=1
F −D
where R = G
. We have,

• F : prediction
• D : data
• G : normalization

It is equivalent to taking uniform prior of X and uncorrelated Gaussian likelihood.

Hence, the algorithm at the posterior pdf is :

X2
log[p(X|D)] = control − (46)
2
Minimizing X 2 is equivalent to minimizing log[p(X|D)]. We take a prior for X and

p(X|D) ∝ p(D|X)p(X) (47)

Example: if we take,
Xj = X0,j + j =⇒ An (48)
uncorrelated Gaussian prior could be used:
M (Xj −X0,j )2
Y 1 22
p(X) = √ e j ∝ e−c/2 (49)
 2π
j=1 j

9
 PM (Xj −X0,j )2
taking c = j=1 2j
.So,

1
L = log(p(X|D)) = constant − (X 2 + c) (50)
2
In this case we need to minimize X 2 + c instead of minimizing X 2 .

12. Describe Linear Partially Observed State Space Models. What is the filter-
ing problem? What is a Bayesian filter?
An Rd stochastic sequence X = {Xk }k≥0 is a state space model if it evolves according
to a recursion of the form:

Xk+1 = Fk Xk + Wk + Uk (51)

Where {Fk }k≥0 is a deterministic sequence of d×d matrices , {Uk }k≥0 is a deterministic
sequence of d × 1 vectors and {Wk }k≥0 is a sequence of independent Rd random vectors
for which EWk = 0 and E[||Wk ||2 ] < ∞.
Partially observed states When are unable to observe directly the state or the
system. However, we observe a corrupted state Xk . In this setting, we assume that
one can directly observe {Zk }k≥0 where:

Zk = Gk Xk + Vk f or k ≥ 0 (52)

Here {Gk }k≥0 is a deterministic sequence of n × d matrices and {Vk }k≥0 is a stochastic
sequence of n × 1 random vectors.
Filtering is concerned with trying to extract the current state Xk from the observed
history of the observation process Z up to time n.

The Bayesian filter gives the complete evolution of the distribution of the state Xk in
terms of the available data in a recursive way. It is done in the following way:
R
1 Time evolution update: π(Xk+1 |{Zj }kj=1 ) = π(Xk+1 |Xk )π(Xk |{Zj }kj=1 ) dXk
π(Zk+1 |Xk+1 )π(Xk+1 |{Zj }kj=1 )
2 Observation update: π(Xk+1 |{Zj }kj=1 ) = π(Zk+1 |{Zj }kj=1 )
R
with: π(Zk+1 |{Zj }kj=1 ) = π(Zk+1 |Xk+1 )π(Xk+1 |{Zj }kj=1 ) dXk+1

13. Describe the Kalman filter, motivate it as a particular case of a Bayesian

filter. Give a simple example of a Kalman Filter applied to a Linear Ob-
served State Space Model.

Given a Bayesian filter for linear evolution and observations subject to Gaussian noise
and initial condition X0 Gaussian, then Xk will be Gaussian. This leads to the Kalman
filter that just tracks the conditional mean and the covariance of Xk .

10
Kalman Filter: Let (Xk , Zk ) follow the linear equation:

Xk+1 = Fk Xk + Wk + Uk (53)

For k ≥ 0 :

(a) {Fk }k≥0 deterministic

(b) {Uk }k≥0
(c) {Wk }k≥0 sequence of independent random vectors Zk = Gk Xk + Vk
(d) {Gk }k≥0 deterministic
(e) {Vk }k≥0 stochastic

We assume that W and V are independent. Then we have,

t
1 Time evolution update: Xk+1|k = Fk+1 X̂k|k + Uk and Σk+1|k = Fk+1 Σk|k Fk+1 +
ΣWk+1 where ΣWk+1 = Cov(Wk+1 )
2 Observation update: X̂k+1|k+1 = X̂k+1|k +Kk+1 (Zk+1 −Gk+1 X̂k+1|k ) and Σk+1|k+1 =
(I − Kk+1 Gk+1 )Σk+1|k
 −1
where: Kk+1 = Σk+1|k Gtk+1 Gk+1Σk+1|k Gtk+1 +ΣVk+1 where ΣVk+1 = Cov(Vk+1 )
such that for n ≥ j :
X̂n|j = E[Xn |Zl : 0 ≤ l ≤ j] (54)
Example: Consider (Ω, F, P ). Given the general stochastic signal evolution for
Un : Ω → R
Z T
−T
Un+1 = Un e + σes−T dWs n = 0, . . . , N − 1 (55)
0

for T = 1 and
1
σ2
Z
σes−1 dWs ∼ N (0, (1 − e−2 )) (56)
0 2
and
Yn+1 = Un+1 + ηn+1 ηn ∼ N (0, T ) (57)
then given an initial condition U0 Gaussian, the filtering distribution,

π(µn+1 |{Yj }n+1

j=1 ) (58)

is Gaussian. So using the Kalman filter we can track the evolution of the first two
moments as follows:

(a) Prediction step:

Un+1 = e−1 m̂n Cn+1 = e−2 Ĉn + Σ

11
 (b) Update step:

m̂n+1 = (I − Kn+1 )mn+1 + Kn+1 Yn+1 Ĉn+1 = (I − Kn+1 )Cn+1

−1
where Kn+1 = Cn+1 Sn+1 and Sn+1 = T + Cn+1
14. Define Markov chains and give an example. Compute its limit distribution
if it exists.

Def. Y = {Yt , t = 0, 1, 2, . . .} is a Markov Chain with state space S, initial distribution

π and transition matrix P if :
(a) Y0 has a distribution π
(b) Conditional distribution of Yt+1 given Yt = i is p(i, i + 1) at time t + 1 and it is
idempotent of Y0 , . . . Yt−1 .
Example:
Every month Raul travels (T ) or not N . Given a state, the transtion probability is
given by ,

p(T → T ) = 0.9 p(T → N ) = 0.1 p(N → T ) = 0.5 p(N → N ) = 0.5 (59)

Then we have the transition matrix,

 
0.9 0.1
P =
0.5 0.5

Then the Markov Chain X ∼ M C(λ, p) is a regular M C such that:

lim P n = (1, . . . , 1)t w̄ (60)

n→∞

where w̄ is the distribution that solves

(
w̄ = w̄P
P (61)
i w̄i = 1

We find the solution to be w̄ = [5/6, 1/6]

15. Let X(t) be a discrete Markov Chain and

u(y, T ) = P (X(T ) = y|X(0) = x).

These probabilities satisfy an evolution equation moving forward in time:

state it and derive it.

Define a µ(x, t) = p(X(t) = x) where x ∈ S sample space. Starting from the numbers
µ(x, 0) = λx ∀x ∈ S and given the transition matrix P . The next evolution in time,

µ(x, t) = λx pt ∀x ∈ S (62)

12
 By conditional probability,

µ(x, t + 1) = p(X(t + 1) = x)
X
= p(x(t + 1) = x|x(t) = y)p(X(t) = y)
y∈S (63)
X
= p(y, x)µ(y, t)
y∈S

We will write the last equation in a matrix form. Suppose that S = {x1 , . . . , xN } and
suppose that |S| < ∞. We define a row vector of order |S|:

µ(t) =< µ1 (x1 , t), . . . , µn (xN , t) > (64)

And so we can write,

u(t + 1) = u(t)P (65)
Where λ = µ(0) =< p(x(0) = x1 ), . . . , p(x(0) = xN ) >

16. Let X(t) be a discrete Markov Chain and for t < T let

f (x, t) = E(V (X(T ))|X(t) = x).

State and derive a backward equation for the above expected value of a
state observable.

Def. F (x, T ) = E[V[X(T )]|X(t) = x] for t < T and x ∈ S. Given F (x, t) =

V (x) ∀x ∈ S and given the transition matrix P , we have the following backward
evolution,
F (x, t) = P T −t f (x, T ) ∀x ∈ S (66)

We derive,

f (x, t) = E[V (X(t))|X(t) = x]

X
= E[V (X|T )|X(t) = x and X(t + 1) = x]p(X(t + 1) = y|X(t) = x)
y∈S (67)
X X
= E[V (X(T ))|X(t + 1) = y]p(x, y) = f (y, t + 1)p(x, y)
y∈S y∈S

The evolution equation of the expectation can be written in a matrix form. Again,
suppose that S = {x1 , . . . , xN } and suppose that |S| < ∞. We define a vector of order
|S|: f (t) =< f1 (t), . . . , fN (t) >. And so we have,

f (t) = P f (t + 1) with f (T ) = V (x) ∀x ∈ S (68)

17. State and discuss a method for sampling from the invariant measure of a
Markov Chain. Mention a possible application of it.
Metropolis-Hasting algorithm samples from pdf f . Fixing f = w , it is possible to

13
sample out of w and then construct an ergodic, recursive MC such that the given w
satisfies , the detailed balance condition.
Metropolis-Hasting:
Input: f the target pdf to be sampled and q(y|x) a proposal density. Let Xn be the
current sample. We perform the following:

1 Given Xn , we sample Yn ∼ q(·|Xn ).

2 Let the acceptance probability be,

f (y)q(x|y)
p(x, y) = min{1, } (69)
f (x)q(y|x)

and take, (
Yn with probability p(Xn , Yn )
Xn+1 = (70)
Xn else

In particular , if q(x|y) = q(y|x) =⇒ p(x, y) = min{1, ff (x)

(y)
}

Next, we check the detailed balance property, we have to show that:

f (x)t(x, y) = f (y)t(y, x) (71)

where t(x, y) are the transition probabilities of the MC defined by Metropolis-Hasting.

The transition pdf is :
Z
t(x, y) = p(x, y)q(y, x) + δx (y)(1 − p(x, s)q(s, x) ds) (72)

By two identities corresponding to two terms in f (x)t(x, y).

(a) ρ(x, y)q(y, x)f (x) = ρ(y, x)q(x|y)f (y).

(b) (1 − r(x))δx (y)f (x) = (1 − r(y))δy (x)f (y)

Hence the detailed balance follows and then f is a stationary pdf for MC.
Application: In the Bayesian framework, we have data {Xi }ni=1 we want to estimate
the parameter θ:
p(θ|{Xi }) ∝ L({Xi }|θ)p(θ) (73)
this has to be known up to a multiplicative constant. Then we can apply M-H algorithm
to draw samples from the posterior p(θ|{Xi }).

14