在scipy中创建新的发行版

Question

我正在尝试根据我拥有的一些数据创建一个分布,然后从该分布中随机绘制.这就是我所拥有的:

from scipy import stats
import numpy

def getDistribution(data):
    kernel = stats.gaussian_kde(data)
    class rv(stats.rv_continuous):
        def _cdf(self, x):
            return kernel.integrate_box_1d(-numpy.Inf, x)
    return rv()

if __name__ == "__main__":
    # pretend this is real data
    data = numpy.concatenate((numpy.random.normal(2,5,100), numpy.random.normal(25,5,100)))
    d = getDistribution(data)

    print d.rvs(size=100) # this usually fails

我觉得这是做什么我也想,但我经常得到一个错误(见下文),当我尝试这样做d.rvs(),并d.rvs(100)永远不会奏效.难道我做错了什么？有更简单或更好的方法吗？如果它是scipy中的一个bug,有没有办法解决它？

最后,是否有更多关于在某处创建自定义发行版的文档？我发现的最好的是scipy.stats.rv_continuous文档,它非常简洁,不包含任何有用的示例.

追溯:

回溯(最近一次调用最后一次):文件"testDistributions.py",第19行,打印d.rvs(size = 100)文件"/usr/local/lib/python2.6/dist-packages/scipy-0.10.0 -py2.6-linux-x86_64.egg/scipy/stats/distributions.py",第696行,在rvs vals = self._rvs(*args)文件"/usr/local/lib/python2.6/dist-packages /scipy-0.10.0-py2.6-linux-x86_64.egg/scipy/stats/distributions.py",第1193行,在_rvs Y = self._ppf(U,*args)文件"/ usr/local/lib /python2.6/dist-packages/scipy-0.10.0-py2.6-linux-x86_64.egg/scipy/stats/distributions.py",第1212行,在_ppf中返回self.vecfunc(q,*args)文件"/usr/local/lib/python2.6/dist-packages/numpy-1.6.1-py2.6-linux-x86_64.egg/numpy/lib/function_base.py",第1862行,致电 theout = self. thefunc(*newargs)文件"/usr/local/lib/python2.6/dist-packages/scipy-0.10.0-py2.6-linux-x86_64.egg/scipy/stats/distributions.py",第1158行,在_ppf_single_call中返回optimize.brentq(self._ppf_to_solve,self.xa,self.xb,args =(q,)+ args,xtol = self.xtol)文件"/usr/local/lib/python2.6/dist-pac kages/scipy-0.10.0-py2.6-linux-x86_64.egg/scipy/optimize/zeros.py",第366行,brentq r = _zeros._brentq(f,a,b,xtol,maxiter,args, full_output,disp)ValueError:f(a)和f(b)必须有不同的符号

编辑

对于那些好奇的人,按照下面答案中的建议,这里的代码有效:

from scipy import stats
import numpy

def getDistribution(data):
    kernel = stats.gaussian_kde(data)
    class rv(stats.rv_continuous):
        def _rvs(self, *x, **y):
            # don't ask me why it's using self._size 
            # nor why I have to cast to int
            return kernel.resample(int(self._size)) 
        def _cdf(self, x):
            return kernel.integrate_box_1d(-numpy.Inf, x)
        def _pdf(self, x):
            return kernel.evaluate(x)
    return rv(name='kdedist', xa=-200, xb=200)

Answer 1

特别是你的追溯:

rvs使用cdf的倒数ppf来创建随机数.由于您没有指定ppf,因此它是通过rootfinding算法计算的brentq.brentq在函数为零时使用下限和上限来搜索值为at(找到x使得cdf(x)= q,q为分位数).

默认的限额,xa并且xb,是在你的例子太少.在创建函数实例时xa,xb可以设置以下适用于我的scipy 0.9.0

def getDistribution(data):
    kernel = stats.gaussian_kde(data)
    class rv(stats.rv_continuous):
        def _cdf(self, x):
            return kernel.integrate_box_1d(-numpy.Inf, x)
    return rv(name='kdedist', xa=-200, xb=200)

目前用于SciPy的改善这种拉请求,所以在接下来的版本xa,并xb会自动扩大,以避免f(a) and f(b) must have different signs例外.

没有太多关于此的文档,最简单的是遵循一些示例(并在邮件列表中询问).

编辑:另外

pdf:由于你也有gaussian_kde给出的密度函数,我会添加_pdf方法,这将使一些计算更有效.

edit2:另外

rvs:如果你有兴趣生成随机数,那么gaussian_kde有一个重采样方法.可以通过从数据中采样并添加高斯噪声来生成随机样本.因此,这将比使用ppf方法的通用rv更快.我会编写一个只调用gaussian_kde的resample方法的._rvs方法.

预计算ppf:我不知道预先计算ppf的任何一般方法.然而,我想这样做的方式(但迄今为止从未尝试过)是在许多点预先计算ppf然后使用线性插值来近似ppf函数.

编辑3:即将_rvs在评论中回答Srivatsan的问题

_rvs是公共方法调用的特定于分发的方法rvs.rvs是一种通用方法,它执行一些参数检查,添加位置和比例,并设置属性self._size,该属性是所请求的随机变量数组的大小,然后调用特定于分发的方法._rvs或它的通用对应方法.额外的参数._rvs是形状参数,但由于在这种情况下没有,*x并且**y是冗余和未使用的.

我不知道size该.rvs方法在多变量情况下的效果如何.这些分布是针对单变量分布而设计的,可能不适用于多变量情况,或者可能需要进行一些重构.