- #1
Graham87
- 71
- 16
- Homework Statement
- Monte Carlo simulation, rejection algorithm
- Relevant Equations
- See pictures
I'm working on this Monte Carlo simulation problem for radiation propagation.
So far I think I have accomplished creating a histogram of simulated values (iii). Now I need to plot it against the exact PDF. I need to find the normalized PDF, but can't obtain the normalization constant, and so on. Or maybe it is working but my histogram is incorrect?
So far I think I have accomplished creating a histogram of simulated values (iii). Now I need to plot it against the exact PDF. I need to find the normalized PDF, but can't obtain the normalization constant, and so on. Or maybe it is working but my histogram is incorrect?
Python:
import numpy as np
import matplotlib.pyplot as plt
def p_kappa_e(kappa_e, kappa):
"""
Computes the value of the unnormalized probability density function p(kappa_c).
Parameters:
kappa_c : float or np.ndarray
The value of kappa_c (can be scalar or array).
kappa : float
The parameter kappa.
Returns:
float or np.ndarray
The value of p(kappa_c).
"""
if np.any(kappa_e >= kappa): # Prevent invalid values for array or scalar input
raise ValueError("kappa_e must be less than kappa for all values.")
term1 = 2
term2 = kappa_e / (kappa * (kappa - kappa_e))
term3 = (kappa_e / (kappa * (kappa - kappa_e))) + (kappa_e - 2)
return term1 + term2 * term3
def kappa_e_max(kappa):
"""
Computes the value of kappa_e,max.
Parameters:
kappa : float
The value of kappa.
Returns:
float
The value of kappa_e,max.
"""
return (2 * kappa**2) / (1 + 2 * kappa)
def p_kappa_e_max(kappa_e_max):
"""
Computes the value of p(kappa_e_max).
Parameters:
kappa_e_max : float
The value of kappa_e,max.
Returns:
float
The value of p(kappa_e_max).
"""
return 2 + 2 * kappa_e_max
# Main parameters
kappa = 0.2
# Parameters for the integral
a_val = 0
b_val = kappa_e_max(kappa)
c = 2
d = p_kappa_e_max(b_val)
# Multiplicative Congruential Generator (MCG)
def mcg_random(seed, a, m, N):
X = np.zeros(N)
X[0] = seed
for i in range(1, N):
X[i] = (a * X[i-1]) % m
return X / m # Normalize to get numbers in the range [0, 1]
# MCG parameters
seed = 12345 # Seed value for random number generator
a = 7**5 # Multiplier
m = 2**31-1 # Modulus
N = 10**6 # Number of points to sample
# Generate random numbers using MCG
x_mcg = mcg_random(seed, a, m, N) * b_val # Scale to [a, b]
y_mcg = mcg_random(seed + 1, a, m, N) * (d - c) + c # Scale to [c, d]
# Rejection sampling: Check how many points lie below the curve p_kappa_c(x)
q_i = y_mcg < p_kappa_e(x_mcg, kappa) # Vectorized operation
# Extract accepted x values
x_accepted = x_mcg[q_i]
# Create the histogram
bin_edges = np.linspace(0, kappa_e_max(kappa), 101) # Define bin edges for 30 bins in (-3, +3)
hist, bins = np.histogram(x_accepted, bins=bin_edges, density=True) # Normalize
# Plot histogram of accepted x values
plt.hist(x_accepted, bins=101, density=True, alpha=0.75, color='blue', edgecolor='black')
plt.title('Histogram of Accepted $\\kappa_c$ Values')
plt.xlabel('$\\kappa_c$')
plt.ylabel('Frequency')
plt.grid(True)
plt.show()
# Calculate mean of q_i
q_mean = np.mean(q_i)