在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文档，它非常简洁并且没有有用的例子。

回溯：

Traceback (most recent call last): File "testDistributions.py", line 19, in print d.rvs(size=100) File "/usr/local/lib/python2.6/dist-packages/scipy-0.10.0-py2.6-linux-x86_64.egg/scipy/stats/distributions.py", line 696, in rvs vals = self._rvs(*args) File "/usr/local/lib/python2.6/dist-packages/scipy-0.10.0-py2.6-linux-x86_64.egg/scipy/stats/distributions.py", line 1193, in _rvs Y = self._ppf(U,*args) File "/usr/local/lib/python2.6/dist-packages/scipy-0.10.0-py2.6-linux-x86_64.egg/scipy/stats/distributions.py", line 1212, in _ppf return self.vecfunc(q,*args) File "/usr/local/lib/python2.6/dist-packages/numpy-1.6.1-py2.6-linux-x86_64.egg/numpy/lib/function_base.py", line 1862, in call theout = self.thefunc(*newargs) File "/usr/local/lib/python2.6/dist-packages/scipy-0.10.0-py2.6-linux-x86_64.egg/scipy/stats/distributions.py", line 1158, in _ppf_single_call return optimize.brentq(self._ppf_to_solve, self.xa, self.xb, args=(q,)+args, xtol=self.xtol) File "/usr/local/lib/python2.6/dist-packages/scipy-0.10.0-py2.6-linux-x86_64.egg/scipy/optimize/zeros.py", line 366, in brentq r = _zeros._brentq(f,a,b,xtol,maxiter,args,full_output,disp) ValueError: f(a) and f(b) must have different signs

编辑

对于那些好奇的，依照下列答案的建议，这里的代码工作：

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，因此它是通过查找算法brentq来计算的。 brentq使用下限和上限，它应该在哪里搜索值，函数为零（找到x使得cdf（x）= q，q是分位数）。

在您的示例中，限制的缺省值xa和xb太小。我与SciPy的0.9.0，xa下面的作品，xb可以在创建函数实例

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方法的通用rvs更快。我会写一个只调用gaussian_kde的resample方法的._rvs方法。

预计算ppf：我不知道任何通用的方法来预先计算ppf。然而，我认为这样做的方式（但从未尝试过）是在多点预先计算ppf，然后使用线性插值来近似ppf函数。

EDIT3：约_rvs回答Srivatsan的问题在评论

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

我不知道在多元情况下size或.rvs方法的形状有多好。这些分布是为单变量分布而设计的，可能不适用于多变量分布情况，或者可能需要一些重构。