ASSIGNMENT: ADVANCED SAMPLING STRATEGIES AND
MARKOV CHAIN FUNDAMENTALS
Physics & Computational Methods
IMPORTANT INSTRUCTIONS
• All questions in this assignment are coding questions.
• SUBMISSION FORMAT:
– Create one Notebook (.ipynb) file
– Solve all questions in this single notebook
– Add clear markdown cells with question titles
– Include comments in your code to explain it wherever needed
– Upload the completed (.ipynb) file to Google Classroom
• For queries contact us on Whatsapp.
Question 1: Monte Carlo Sampling of a Localized Signal
Physical System:
A detector measures signals originating from a localized source centered at the origin. Due to instrumental
effects, the spatial intensity profile of the signal is well approximated by
2
/(2σ 2 )
I(x) = e−x
over the detector range x ∈ [−L, L], where σ characterizes the width of the signal.
The total detected signal strength is defined as
Z L
S= I(x) dx
−L
This integral must be evaluated numerically.
Part (a): Accept-Reject Sampling
Interpret I(x) as an unnormalized probability density.
You may take L = 3.
Use a uniform proposal distribution
1
g(x) = for x ∈ [−L, L]
2L
• Determine the envelope constant M such that
I(x) ≤ M g(x) for all x ∈ [−L, L]
1
2
• Implement accept-reject sampling to generate N = 50,000 samples distributed according to I(x).
• Store the accepted samples and compute the acceptance rate.
• Plot a normalized histogram of the accepted samples and overlay the function I(x).
Part (b): Estimating the Total Signal
The total signal strength is
Z L
S= I(x) dx
−L
Show that the integral can be estimated as
S ≈ (2L)M × (acceptance rate)
where the acceptance rate is the ratio of number of accepted samples to the value of N.
Explain briefly why the acceptance rate equals the fraction of the rectangular area under the curve I(x).
(Hint: Think graphically if needed.)
Part (c): Signal Width
Define the mean squared position
Z L
1
⟨x2 ⟩ = x2 I(x) dx
S −L
2
Estimate ⟨x ⟩ using the Monte Carlo samples. (Compute the average using the samples gathered.)
Physical Interpretation: Explain in 2 to 3 sentences how the choice of proposal distribution affects
sampling efficiency and how accept-reject sampling avoids explicit normalization.
Question 2: Markov Chain Relaxation and Stationary Distributions
Physical System:
Consider a physical system that can exist in three discrete states due to thermal fluctuations. Each state
corresponds to a different energy level:
• State A: EA = 0 (ground state)
• State B: EB = 5 (intermediate state)
• State C: EC = 12 (excited state)
At a fixed temperature, the system evolves stochastically between these states. The dynamics are modeled
as a discrete-time Markov chain with transition matrix
0.85 0.10 0.05
P = 0.15 0.75 0.10
0.30 0.20 0.50
where Pij denotes the probability of transitioning from state i to state j in one timestep. The diagonal entries
represent the probability of remaining in the same state. These are governed by Boltzmann distributions.
The probability distribution over states evolves according to
π t+1 = π t · P
3
Part (a): Stationary Distribution
A stationary distribution π stat satisfies
π stat · P = π stat
Construct the matrix P in Python and compute the stationary distribution using [Link](P.T).
Identify the eigenvector corresponding to eigenvalue 1 and normalize it so that the probabilities sum to 1.
Print the stationary probabilities in the format:
A = [Link], B = [Link], C = [Link]
Physical Interpretation: Explain in 2 to 3 sentences why the lowest-energy state does not necessarily
have probability 1, and how the transition probabilities influence the stationary distribution. (Refer to solved
example if needed.)
Part (b): Time Evolution and Thermal Relaxation
Assume the system starts entirely in the ground state:
π 0 = [1, 0, 0]
Simulate the time evolution of the system for 500 timesteps by iterating
π t+1 = π t · P
using a loop and the @ operator for matrix multiplication.
Store the probabilities of each state at every timestep and plot them as functions of time. Add horizontal
dashed lines at the stationary values from part (a). Include proper axis labels and a legend.
Interpretation: In 2 to 3 sentences, describe how quickly the system reaches equilibrium and which
state equilibrates first. Relate your answer to the structure of the transition matrix.
Question 3: Importance Sampling for Spectral Power Estimation
Physical System:
Consider a thermal radiation source that emits power across a spectrum of frequencies. The spectral
power density (power per unit frequency) is given by:
f2
P (f ) =
ef /T −1
where f ∈ [0, 50] GHz is the frequency and T = 10 is a temperature-like parameter (in normalized units).
This is related to the Planck distribution for blackbody radiation.
The total radiated power is: Z 50
Φ= P (f ) df
0
This integral cannot be evaluated analytically and must be estimated numerically using Monte Carlo
methods.
Part (a): Inspect the Integrand and Choose a Proposal
Examine the integrand by plotting it:
f2
P (f ) =
ef /T −1
4
• What is its behavior at small frequencies? (Hint: Does it grow, decay, or vanish?)
• What is its behavior at large frequencies? (Hint: How does the denominator dominate?)
• Where is the peak? (Estimate numerically or graphically.)
Based on your observations, choose a proposal distribution q(f ) that:
1. Peaks near the maximum of P (f )
2. Has a similar decay at large frequencies
3. Is easy to sample from (e.g., uniform, triangular, or exponential on a truncated domain)
Justify your choice: In 3-4 sentences, explain why your proposal distribution is a sensible approxi-
mation to the envelope of P (f ), and how it will reduce variance compared to naive uniform sampling over
[0, 50].
Part (b): Implement Two Monte Carlo Estimates
1. Naive Estimate: Sample N = 50,000 frequencies uniformly from [0, 50]:
1
qnaive (f ) =
50
Estimate the integral as:
N
1 X
Φnaive ≈ P (fi ) × 50
N i=1
Record the estimate and compute the variance of the samples.
2. Importance Sampling Estimate: Sample N = 50,000 frequencies from your chosen proposal q(f ).
Estimate:
N
1 X P (fi )
ΦIS ≈
N i=1 q(fi )
Record the estimate and compute the variance of the samples.
The variance of each Monte Carlo estimate is computed as:
N
1 X 2
Var = Xi − X
N i=1
where Xi are the individual weighted samples and X is their mean. Alternatively, use [Link]() on the
array of weighted samples.
Part (c): Variance Reduction and Physical Insight
Physical Interpretation: In 3-4 sentences, explain:
• Why does the naive method waste samples in regions where P (f ) is negligible?
• How does matching the proposal to the envelope of P (f ) concentrate samples where they matter? Why
did you choose that specific proposal?
5
HINTS
• Q1: (See Question 1) For part (a), use [Link]()
h −x2 /(2σ2 ) i to create a fine grid over [−L, L] and com-
e
pute I(x) at each point. Find M as maxx 1/(2L) . For accept-reject, generate uniform proposals
with [Link](), compute acceptance probabilities, and use [Link]() to decide
acceptance. For the histogram, use [Link]() and normalize by dividing by N × ∆x.
• Q2: (See Question 2) For part (a), use [Link](P.T) to find eigenvalues and eigenvectors.
The eigenvector corresponding to eigenvalue 1 (or closest to 1) is the stationary distribution. Normalize
it by dividing by its sum. For part (b), use a loop with the @ operator: pi t = pi t @ P. Store values
in arrays before updating pi t.
• Q3: (See Question 3) For part (a), plot P (f ) over [0, 50] to visually identify the peak and behavior.
For part (b), a good proposal might be exponential decay or triangular distribution peaking near
f ≈ 25. Use inverse CDF sampling for non-uniform proposals. For variance calculation, create arrays
of weighted samples and use [Link]().