单能伽马射线模拟：交互深度采样及平均自由路径代码正确性验证

嘿，我一眼就看出你当前的代码存在几个关键问题，完全没模拟出光电吸收的交互深度分布，下面给你拆解并修正：

核心问题所在 #

• 错误的分布选择：你用了np.random.randn()生成正态分布随机数，但光电吸收的交互深度应该服从指数分布——这是因为伽马光子在介质中传播时，未发生交互就到达深度x的概率是 ( P(x) = e^{-\mu x} )，对应的采样分布就是指数分布。
• 随机种子设置错误：random.seed=()这种写法不合法，正确的种子设置应该是np.random.seed(42)（42是常用的固定种子，方便复现结果）或者不填参数用系统随机种子。
• 数据生成逻辑完全偏离需求：你当前生成的data = np.random.randn(50000)*10是均值为0、标准差为10的正态分布数据，和你需要的平均自由路径0.5cm（( 1/\mu = 1/2 = 0.5 )）完全不匹配。

修正后的代码 #

import matplotlib.pyplot as plt
import numpy as np

# 设置参数
mu = 2  # 线性衰减系数，单位cm⁻¹
mean_free_path = 1 / mu  # 平均自由路径，计算得0.5cm
sample_size = 50000

# 设置随机种子，保证结果可复现
np.random.seed(42)

# 生成指数分布的交互深度样本
# np.random.exponential的scale参数对应分布的均值，也就是我们的平均自由路径
interaction_depths = np.random.exponential(scale=mean_free_path, size=sample_size)

# 绘制直方图，同时添加理论概率密度曲线
bins = np.arange(0, interaction_depths.max() + 0.1, 0.1)
plt.hist(interaction_depths, bins=bins, density=True, alpha=0.7, label='模拟样本')

# 生成理论概率密度曲线数据
x = np.linspace(0, interaction_depths.max(), 1000)
pdf = mu * np.exp(-mu * x)
plt.plot(x, pdf, 'r-', label='理论PDF')

# 图表美化
plt.xlabel('交互深度 (cm)')
plt.ylabel('概率密度')
plt.title('光电吸收交互深度分布（模拟 vs 理论）')
plt.legend()
plt.show()

# 输出模拟样本的均值，和理论值对比
print(f"理论平均自由路径: {mean_free_path:.2f} cm")
print(f"模拟样本均值: {np.mean(interaction_depths):.2f} cm")

验证逻辑说明 #

运行修正后的代码，你会看到：

• 模拟样本的直方图和红色的理论概率密度曲线几乎重合
• 模拟样本的均值会非常接近0.5cm（随着样本量增大，误差会越来越小），这就验证了模拟的正确性

内容的提问来源于stack exchange，提问作者James