カスタムPDF scipy.stats.rv_continuous不要な上限

Question

から、私はQSOの形で一定の明るさのランダムな確率密度関数を生成しようとしています：

1 /（（L/L_B^*）^アルファ（L/L_B^*）^β）

ここで、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のパワーにrasedれる）

分布を生成することになっていますthisに近づき、無限遠まで下がるグラフですが、実際にはグラフはthis（10,000サンプル）のようになります。上限は、使用されるサンプルの量に関係なく同じです。制限されている理由はありますか？

Answer 1

あなたのPDFは正しく正規化されていません。

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の積分は1.あなたのPDFは約3.4712に統合しなければなりません。それだけで、あなたのケースで、実際には、あなたがヒストグラムに表示されているものについて2.325

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

である、0からxにPDFの積分が最初に1.0に達したところまでのドメインからサンプルを生成します。 。

は、問題を検証するために簡単な修正として、私はに_pdf()方法のreturnステートメントを変更：

 return D/3.47121849657489 

て、もう一度スクリプトを実行しました。 （実際の修正では、その値は他のパラメータの関数となります。）次にコマンド

In [85]: import matplotlib.pyplot as plt 

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

はこのプロットを生成します。