Monte Carlo Simulation（Rejection Algorithm）求助：归一化PDF求解与直方图修正

我是一名物理系学生，正在处理一个蒙特卡洛模拟的问题。（相关问题示意图及归一化常数公式图已省略）

目前我应该已经完成了模拟值的直方图绘制，但接下来需要将其与精确的PDF进行对比。我需要得到归一化后的PDF，却卡在了归一化常数的求解上；另外我的直方图似乎也存在问题，尝试修正但还没找到头绪。

以下是我目前写的代码：

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)

备注：内容来源于stack exchange，提问作者Graham_PP