2

次のような Schechter Luminosity 関数を扱っています。

phi(L)dL = norm. Factor * (L/Lstar)^(a) * exp (L/Lstar) d(L/Lstar)

たとえば、L/Lstar は l です。

の累積分布関数の解析解は、ガンマ関数によって与えられます: N = ノルム係数*ガンマ(a+1, l)。

積分の限界が L から無限大であるため、これは不完全なガンマ関数です。

現在、Python で cdf をプロットしようとしています。私が使用した：

import scipy.special as ss
si= [ss.gammainc(a+1, l[i]) for i in a]  #cdf

(ここで、l[I] は乱数から作成した配列です)

結果のグラフの合計は 1 で、cdf のように見えます。しかし今、私はそれをランダム化したいです。したがって、cdf = 1 の代わりに、cdf = 乱数を設定します (Python によって一様に生成されます)。ここで、ランダム サンプリングを使用してカウント数と L のヒストグラムをプロットする場合は、ガンマ関数を反転する必要があります。

私の質問は: Python でガンマ関数を反転するにはどうすればよいですか?

これは私が今持っているものです：

u= [random.uniform(0,1) for i in a]

l= [ss.gammaincinv(a+1, u[i]) for i in a]

plt.plot(l, u, '.')

plt.show()

plt.hist(l, bins=50,rwidth= 1.5,histtype='step', lw= 0.7, normed= True, range=(-0.5, 1))

plt.show()

コンパイラは文句を言いませんが、ヒストグラムは間違った形です。cdf のランダムにサンプリングされたヒストグラムは、PDF の形状を回復するはずだと思います。

私は何を間違っていますか？どうやら、不完全なガンマ関数の scipy のバージョンは「正則化」されています。つまり、完全なガンマ関数で除算されます。したがって、 gammainc(a+1, u[I])* gamma(a+1) を掛けても、まだ機能しません。

軸は対数スケールです。

助言がありますか？

結論: ランダム サンプリングにより、Schechter 光度関数の cdf のヒストグラムを作成する必要があります。

4

1 に答える 1

1

初挑戦：

関数は、ドメインから範囲へのマッピングです。したがって、次のように記述できます。

def function(x):
    # ...

Domain = list(range(0, 1000)) # [0,1000)
mapping = {}
inverse_mapping = {}
for x in Domain:
    y = function(x)
    mapping[x] = y
    inverse_mapping[y] = x

def inverse_function(y):
    return inverse_mapping[y] # not a continuous function. needs improvement

このようなことが頭に浮かんだら、私に知らせてください。cdf のような単調な関数については改善できます。

于 2013-11-20T20:19:57.927 に答える