python - SciPy:从 PMF 生成自定义随机变量

我试图在 Python 中根据某个丑陋的分布生成随机变量。我有一个 PMF 的显式表达式，但它涉及一些产品，这使得获取和反转 CDF 令人不快(请参阅下面的代码以了解 PMF 的显式形式)。

本质上，我试图通过其 PMF 在 Python 中定义一个随机变量，然后让内置代码完成从分布中采样的艰苦工作。如果 RV 的支持是有限的，我知道如何做到这一点，但这里的支持是可数无限的。

我目前尝试按照以下@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

基本上，似乎 np.arange() 列表在 def _pmf() 函数内部不起作用。我不知道为什么。任何人都可以在这里启发我和/或指出解决方法吗？

编辑 1: 解决了 askewchan 的一些问题，上面反射(reflect)了编辑。

编辑 2: askewchan 使用阶乘函数提出了一个有趣的近似值，但我正在寻找更多的精确解决方案，例如我正在尝试使用 np.arange 的解决方案。

最佳答案

您应该能够像这样子类化 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)

关于python - SciPy:从 PMF 生成自定义随机变量，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/23276631/