python - Python 中具有 Gamma 函数的累积分布函数

Question

我正在处理如下所示的 Schechter Luminosity 函数:

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

比方说，L/Lstar 是 l。

其累积分布函数的解析解由 Gamma 函数给出: N = 范数因子* Gamma(a+1, l)。

这是不完全 Gamma 函数，因为积分的极限是 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 = 随机数(由 Python 统一生成)，而不是 cdf = 1。现在，如果我想通过随机抽样绘制计数与 L 的直方图，我需要反转 Gamma 函数。

我的问题是:如何在 Python 中反转 Gamma 函数？

这是我现在拥有的:

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 版本的不完整 Gamma 函数是“正则化”的，这意味着它被完整的 Gamma 函数除以。所以如果我乘以 gammainc(a+1, u[I])* gamma(a+1) 它仍然不起作用。

轴是对数刻度的。

有什么建议吗？

底线:我需要通过随机抽样制作 Schechter 光度函数的 cdf 的直方图。

Answer 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 这样的单调函数对其进行改进。