我正在处理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的直方图 .