Моделирование Монте-Карло, алгоритм отклонения ⇐ Python

[Anonymous]

Я студент-физик, работаю над задачей моделирования Монте-Карло.
Задача моделирования Монте-Карло
Формула константы нормализации
На данный момент я думаю, что мне удалось создать гистограмму смоделированных значений. Теперь мне нужно сопоставить его с точным 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 = (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)


Подробнее здесь: https://stackoverflow.com/questions/793 ... -algorithm