python - 在 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 中的错误，有什么办法可以解决它吗？

最后，是否有更多关于在某处创建自定义发行版的文档？我找到的最好的是 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 的 pull request 来改进这个，所以在下一个版本中 xa 和 xb 将自动扩展以避免 f(a ) 和 f(b) 必须有不同的符号异常(exception)。

这方面的文档不多，最简单的方法是遵循一些示例(并在邮件列表中询问)。

编辑:添加

pdf:既然你有gaussian_kde也给出的密度函数，我会添加_pdf方法，这将使一些计算更有效率。

edit2: 添加

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

预计算 ppf:我不知道有什么通用的方法可以预计算 ppf。然而，我想到的方法(但到目前为止从未尝试过)是在许多点预先计算 ppf，然后使用线性插值来近似 ppf 函数。

edit3: 关于 _rvs 回答 Srivatsan 在评论中的问题

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

我不知道 size 或 .rvs 方法的形状在多变量情况下的效果如何。这些分布是为单变量分布设计的，可能无法完全适用于多变量情况，或者可能需要进行一些整形。