User: How to use Metropolis Algorithm to do Monte Carlo summation for a integral : int x^2 exp(-x^2/2) dx?
Assistant: To use the Metropolis algorithm for Monte Carlo integration of the function \( f(x) = x^2 \exp\left(-\frac{x^2}{2}\right) \), you can follow these steps:
We are trying to estimate the integral:
\[ I = \int_{-\infty}^{\infty} x^2 \exp\left(-\frac{x^2}{2}\right) dx \]
This resembles a Gaussian integral, but Monte Carlo methods are useful when direct evaluation is difficult or when sampling is preferred.
The Metropolis algorithm generates samples according to a target distribution (often a probability density function). We can use it to generate samples from the distribution and approximate the integral as an average over those samples.
The integral looks similar to a Gaussian distribution with a variance of 1. So we can interpret the exponential term as part of a normal distribution:
\[ p(x) = \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{x^2}{2}\right) \]
Thus, we want to sample from this distribution and compute the expectation value of \(x^2\), as:
\[ I = \sqrt{2\pi} \cdot \mathbb{E}_{p(x)}[x^2] \]
\[ A(x' | x) = \min \left( 1, \frac{p(x')}{p(x)} \right) \]
Since \(p(x) = \exp\left(-\frac{x^2}{2}\right)\), the acceptance ratio simplifies to:
\[ A(x' | x) = \min \left( 1, \exp\left(-\frac{x'^2 - x^2}{2}\right) \right) \]
\[ \mathbb{E}[x^2] \approx \frac{1}{N} \sum_{i=1}^N x_i^2 \]
\[ I \approx \sqrt{2\pi} \cdot \frac{1}{N} \sum_{i=1}^N x_i^2 \]
import numpy as np
# Target distribution: exp(-x^2 / 2)
def target_distribution(x):
return np.exp(-x**2 / 2)
# Metropolis algorithm
def metropolis_hastings(N, sigma):
x = 0 # initial guess
samples = []
for _ in range(N):
# Propose new state from a Gaussian around current state
x_new = x + np.random.normal(0, sigma)
# Acceptance ratio
acceptance_ratio = min(1, target_distribution(x_new) / target_distribution(x))
# Accept or reject
if np.random.rand() < acceptance_ratio:
x = x_new
samples.append(x)
return np.array(samples)
# Parameters
N = 100000 # Number of samples
sigma = 1.0 # Step size in the proposal distribution
# Generate samples
samples = metropolis_hastings(N, sigma)
# Compute the expectation value of x^2
expectation_x2 = np.mean(samples**2)
# Multiply by the normalization factor sqrt(2*pi)
integral_estimate = np.sqrt(2 * np.pi) * expectation_x2
print("Estimated value of the integral:", integral_estimate)
User: Please include the exact integral value in the Python code.
Assistant: You can compute the exact value of the integral \( I = \int_{-\infty}^{\infty} x^2 \exp\left(-\frac{x^2}{2}\right) dx \) analytically. This is a well-known Gaussian integral and can be solved using integration by parts or a standard result for Gaussian integrals. The exact value is:
\[ I_{\text{exact}} = \sqrt{2\pi} \]
Here is the updated Python code:
import numpy as np
# Target distribution: exp(-x^2 / 2)
def target_distribution(x):
return np.exp(-x**2 / 2)
# Metropolis algorithm
def metropolis_hastings(N, sigma):
x = 0 # initial guess
samples = []
for _ in range(N):
# Propose new state from a Gaussian around current state
x_new = x + np.random.normal(0, sigma)
# Acceptance ratio
acceptance_ratio = min(1, target_distribution(x_new) / target_distribution(x))
# Accept or reject
if np.random.rand() < acceptance_ratio:
x = x_new
samples.append(x)
return np.array(samples)
# Parameters
N = 100000 # Number of samples
sigma = 1.0 # Step size in the proposal distribution
# Generate samples
samples = metropolis_hastings(N, sigma)
# Compute the expectation value of x^2
expectation_x2 = np.mean(samples**2)
# Multiply by the normalization factor sqrt(2*pi)
integral_estimate = np.sqrt(2 * np.pi) * expectation_x2
# Exact integral value
integral_exact = np.sqrt(2 * np.pi)
# Display results
print(f"Estimated value of the integral: {integral_estimate}")
print(f"Exact value of the integral: {integral_exact}")
print(f"Relative error: {(integral_estimate - integral_exact) / integral_exact:.5f}")
The output would look something like this after running the code:
Estimated value of the integral: 2.50743
Exact value of the integral: 2.50663
Relative error: 0.00032
User: Can you change the way to use choose a new x? Instead of random normal, can we use flat random number?
Assistant: Yes! Instead of proposing a new \( x' \) using a Gaussian distribution (normal distribution), you can propose the new \( x' \) using a uniform distribution over a certain range. The rest of the Metropolis-Hastings algorithm remains the same, but the proposal distribution will now be uniform.
Here is the updated Python code that uses a uniform proposal distribution:
import numpy as np
# Target distribution: exp(-x^2 / 2)
def target_distribution(x):
return np.exp(-x**2 / 2)
# Metropolis algorithm with uniform proposal distribution
def metropolis_hastings_uniform(N, delta):
x = 0 # initial guess
samples = []
for _ in range(N):
# Propose new state from a uniform distribution around current state
x_new = x + np.random.uniform(-delta, delta)
# Acceptance ratio
acceptance_ratio = min(1, target_distribution(x_new) / target_distribution(x))
# Accept or reject
if np.random.rand() < acceptance_ratio:
x = x_new
samples.append(x)
return np.array(samples)
# Parameters
N = 100000 # Number of samples
delta = 2.0 # Step size in the proposal distribution (uniform width)
# Generate samples
samples = metropolis_hastings_uniform(N, delta)
# Compute the expectation value of x^2
expectation_x2 = np.mean(samples**2)
# Multiply by the normalization factor sqrt(2*pi)
integral_estimate = np.sqrt(2 * np.pi) * expectation_x2
# Exact integral value
integral_exact = np.sqrt(2 * np.pi)
# Display results
print(f"Estimated value of the integral: {integral_estimate}")
print(f"Exact value of the integral: {integral_exact}")
print(f"Relative error: {(integral_estimate - integral_exact) / integral_exact:.5f}")
After running the code with a uniform proposal, you should see output similar to this:
Estimated value of the integral: 2.50972
Exact value of the integral: 2.50663
Relative error: 0.00123