Monte Carlo simulation, rejection algorithm

[Graham87]

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?

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)

 

[FactChecker]

Note 2 hints that the normalizing constant should be $p_{MC}(0)/2$.

 

[Graham87]

Note 2 hints that the normalizing constant should be $p_{MC}(0)/2$.

Yes, I know that. I am struggling to code for that. Not sure what $p_{MC}(0)/2$, since it’s not an equation, but a simulated histogram of counts.
And is my code correct so far?

Update: I realised there is something wrong with my histogram, and currently working on it.