SciPy：从PMF生成自定义随机变量

Question

我正在尝试在Python中生成一些随机变量，这些变量遵循一种比较复杂的分布。我有一个明确的概率质量函数（PMF）的表达式，但里面有一些乘法运算，这让求出和反转累积分布函数（CDF）变得很麻烦（下面的代码展示了PMF的具体形式）。

简单来说，我想在Python中通过PMF来定义一个随机变量，然后让内置的代码来处理从这个分布中抽样的工作。如果随机变量的取值是有限的，我知道怎么做，但在这里，取值是可数无限的。

我目前根据@askewchan的建议尝试运行的代码是：

import scipy as sp
import numpy as np

class x_gen(sp.stats.rv_discrete):
    def _pmf(self,k,param):
        num = np.arange(1+param, k+param, 1)
        denom = np.arange(3+2*param, k+3+2*param, 1)

        p = (2+param)*(np.prod(num)/np.prod(denom))

        return p

pa_limit = limitrv_gen()
print pa_limit.rvs(alpha,n=1)

然而，运行时出现了错误：

File "limiting_sim.py", line 42, in _pmf
    num = np.arange(1+param, k+param, 1)
TypeError: only length-1 arrays can be converted to Python scalars

基本上，看起来在def _pmf()函数内部，np.arange()这个列表似乎没有正常工作。我不知道为什么会这样。有没有人能帮我解答一下，或者指出一个解决办法？

编辑 1：根据askewchan的提问，澄清了一些问题，以上的编辑已反映出来。

编辑 2：askewchan建议使用阶乘函数进行一种有趣的近似，但我更想要一个精确的解决方案，比如我正在尝试用np.arange实现的那种。

Answer 1

你可以像下面这样来扩展 rv_discrete 类：

class mydist_gen(rv_discrete):
    def _pmf(self, n, param):
        return yourpmf(n, param)

然后你可以用下面的方式创建一个分布实例：

mydist = mydist_gen()

接着可以用下面的方式生成样本：

mydist.rvs(param, size=1000)

或者你也可以创建一个固定的分布对象，方法是：

mydistp = mydist(param)

最后可以用下面的方式生成样本：

mydistp.rvs(1000)

---

根据你的例子，这样应该是可以工作的，因为 factorial 会自动处理不同的输入。但如果 alpha 的值足够大，可能会出现问题：

import scipy as sp
import numpy as np
from scipy.misc import factorial

class limitrv_gen(sp.stats.rv_discrete):
    def _pmf(self, k, alpha):
        #num = np.prod(np.arange(1+alpha, k+alpha))
        num = factorial(k+alpha-1) / factorial(alpha)
        #denom = np.prod(np.arange(3+2*alpha, k+3+2*alpha))
        denom = factorial(k + 2 + 2*alpha) / factorial(2 + 2*alpha)

        return (2+alpha) * num / denom

pa_limit = limitrv_gen()
alpha = 100
pa_limit.rvs(alpha, size=10)