来自 scipy.stats.rv_continuous 的自定义 PDF 不需要的上限答案

【问题标题】：Custom PDF from scipy.stats.rv_continuous unwanted upper-bound来自 scipy.stats.rv_continuous 的自定义 PDF 不需要的上限
【发布时间】：2017-07-25 15:45:03
【问题描述】：

我正在尝试使用以下形式生成具有一定亮度的 QSO 的随机概率密度函数：

1/( (L/L_B^* )^alpha + (L/L_B^* )^beta )

其中 L_B^*、alpha 和 beta 都是常数。为此，使用以下代码：

import scipy.stats as st

logLbreak = 43.88
alpha = 3.4
beta = 1.6


class my_pdf(st.rv_continuous):

    def _pdf(self,l_L): 
        #"l_L" in this is always log L        
        L = 10**(l_L/logLbreak)
        D = 1/(L**alpha + L**beta)
        return D

dist_Log_L = my_pdf(momtype = 0, a = 0,name='l_L_dist')


distro = dist_Log_L.rvs(size = 10000)

（L/L^* 被提升到 10 次方，因为一切都是在对数刻度中完成的）

该分布应该生成一个近似于this 的图形，逐渐减小到无穷大，但实际上它生成的图形看起来像this（10,000 个样本）。无论使用多少样本，上限都是相同的。是否有理由限制它的方式？

【问题讨论】：

标签： python pdf scipy statistics astronomy

【解决方案1】：

您的 PDF 未正确规范化。 PDF 在域上的积分必须为 1。您的 PDF 积分约为 3.4712：

In [72]: from scipy.integrate import quad

In [73]: quad(dist_Log_L._pdf, 0, 100)
Out[73]: (3.4712183965415373, 2.0134487716044682e-11)

In [74]: quad(dist_Log_L._pdf, 0, 800)
Out[74]: (3.4712184965748905, 2.013626296581202e-11)

In [75]: quad(dist_Log_L._pdf, 0, 1000)
Out[75]: (3.47121849657489, 8.412130378805368e-10)

这将破坏该类对inverse transform sampling 的实现。它只会从域中生成样本，直到 PDF 从 0 到 x 的积分首先达到 1.0，在您的情况下约为 2.325

In [81]: quad(dist_Log_L._pdf, 0, 2.325)
Out[81]: (1.0000875374350238, 1.1103202107010366e-14)

实际上，就是您在直方图中看到的内容。

作为验证问题的快速修复，我将_pdf() 方法的return 语句修改为：

        return D/3.47121849657489

然后再次运行您的脚本。（在真正的修复中，该值将是其他参数的函数。）然后是命令

In [85]: import matplotlib.pyplot as plt

In [86]: plt.hist(distro, bins=31)

生成此图：

【讨论】：

非常感谢！我很确定我使用的常量不正确，所以我想如果我真的正确实现它们，pdf 将是一个真正的 pdf。