定制PDF scipy.stats.rv_continuous不需要的上限

我试图产生某些光度与以下形式的QSO的随机概率密度函数：定制PDF scipy.stats.rv_continuous不需要的上限

1 /（（L/L_B^*）^阿尔法+（L/L_B^*）^β）

其中L_B^*，α和β都是常数。要做到这一点，下面的代码被用于：

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)

（^ *是rased以10的倍数，因为一切都在数比例正在做L/L）

的分布应该产生一个近似于this的图表，拖尾到无穷大，但实际上它生成的图形看起来像this（10,000个样本）。无论使用的样本数量如何，上限都是相同的。是否有理由限制它的方式？

来源

2017-07-25 Matteo Mcdonnell

您的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)

产生这样的情节：

来源

2017-07-25 17:37:24

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

定制PDF scipy.stats.rv_continuous不需要的上限

回答

相关问题