定製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不需要的上限

回答

相關問題