如何绘制Gamma(1,1)分布变量变换后z=x²/(1+sin(x))的概率密度函数

实现方案 #

你要计算非线性变换后变量z的概率密度，最简便的实现方式是蒙特卡洛采样+核密度估计，因为你的变换z=x²/(1+sin(x))不是单调函数，无法直接用一元随机变量变换的PDF公式计算，采样法不需要考虑单调性问题，实现门槛更低。

你原有代码的问题 #

• 你用np.linspace生成的是等间隔x值，不是Gamma分布的样本，不能直接用来推导变换后的z的分布
• stats.gamma.pdf(x, a=1, loc=1)里的loc=1参数错误，Gamma(1,1)分布不需要平移，loc应设为默认值0
• 标签里的α=29, β=3和你需要的Gamma(1,1)分布参数不匹配，属于错误标注

完整可运行代码 #

import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt

# 1. 生成大量服从Gamma(1,1)的随机样本，样本量越大PDF拟合精度越高
x_samples = stats.gamma.rvs(a=1, scale=1, size=100000)
# 2. 按变换公式计算z，加极小值避免sin(x)=-1时分母为0
z_samples = x_samples ** 2 / (1 + np.sin(x_samples) + 1e-8)
# 可选：过滤极端大的z值，避免绘图时显示异常
z_samples = z_samples[z_samples < 500]

# 3. 绘制z的PDF
plt.figure(figsize=(10, 6))
# 绘制归一化直方图作为参考
plt.hist(z_samples, bins=150, density=True, alpha=0.3, color='skyblue', label='采样分布直方图')
# 用高斯核密度估计得到平滑的PDF曲线
kde = stats.gaussian_kde(z_samples)
z_grid = np.linspace(z_samples.min(), z_samples.max(), 500)
plt.plot(z_grid, kde(z_grid), 'r-', linewidth=2, label='KDE拟合PDF')

plt.xlabel('z')
plt.ylabel('概率密度')
plt.xlim(0, 200) # 可根据需要调整x轴显示范围
plt.legend()
plt.show()

其他说明 #

如果需要更精确的理论PDF，可以对变换的分段单调区间分别用变量替换公式求和，但实现复杂度很高，普通场景下采样法的精度完全够用。

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