首页 文章

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

提问于
浏览
3

我试图在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() 函数内部无法正常工作 . 我不知道为什么 . 任何人都可以在这里启发我和/或指出修复?

EDIT 1: 清除了askewchan的一些问题,上面反映的编辑 .

EDIT 2: askewchan建议使用阶乘函数进行有趣的近似,但是我试图使用np.arange .

1 回答

  • 2

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

相关问题