2

我正在处理看起来像这样的 Schechter Luminosity 函数：

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

说，L/Lstar 是 l。

其累积分布函数的解析解由 gamma 函数给出：N = norm factor* 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 的不完全伽马函数版本是“正则化”的，这意味着它被完整的伽马函数划分。因此，如果我将 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 回答