首页 文章

Python中使用Gamma函数的累积分布函数

提问于
浏览 次
2

我正在处理Schechter Luminosity功能，如下所示：

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

说，L / Lstar是l .

其累积分布函数的解析解由伽马函数给出：N =范数因子* Gamma（a 1,1） .

这是不完整的伽马函数，因为积分的极限是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的直方图，随机采样，我需要反转伽玛函数 .

我的问题是： How do I invert the Gamma function in 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（1，u [I]）* gamma（a 1）相乘，它仍然不起作用 .

轴是按比例缩放的 .

有什么建议？

结论：我需要通过随机抽样制作Schechter光度函数cdf的直方图 .

python gamma-function

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这样的单调函数 .

    回复于

相关问题