Monte Carlo simulation, rejection algorithm

Question

I'm a physics student working on a Monte Carlo simulation problem. Monte Carlo Simulation Problem

So far I think I have accomplished creating a histogram of simulated values. 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.

There is also something wrong with my histogram. Trying to fix it but not sure how.

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)

Monte Carlo simulation, rejection algorithm

Answers (1)

Related Questions