在 scipy 中创建新的发行版

我试图根据我拥有的一些数据创建一个分布，然后从该分布中随机抽取。这是我所拥有的：

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(), and d.rvs(100)从来不工作。难道我做错了什么？有没有更简单或更好的方法来做到这一点？如果这是 scipy 中的错误，有什么方法可以解决它吗？

最后，是否有更多关于创建自定义发行版的文档？我发现的最好的文档是 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路，房车 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 行，在calltheout = 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-packages/scipy-0.10.0-py2.6-linux-x86_64.egg/scipy/optimize/zeros.py”， 366 号线，布伦特克 r = _zeros._brentq(f,a,b,xtol,maxiter,args,full_output,disp) ValueError: f(a) 和 f(b) 必须具有不同的符号

Edit

对于那些好奇的人，请按照下面答案中的建议，以下是有效的代码：

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)

特别是你的回溯：

rvs 使用 cdf 的逆 ppf 来创建随机数。由于您没有指定 ppf，因此它是通过求根算法计算的，brentq. brentq使用下限和上限来搜索函数为零的值 at 的位置（找到 x 使得 cdf(x)=q，q 是分位数）。

默认的限制，xa and 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 and xb将自动扩展以避免f(a) and f(b) must have different signs例外。

关于此的文档并不多，最简单的是遵循一些示例（并在邮件列表上询问）。

编辑：添加

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

编辑2：添加

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

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

编辑3：关于_rvs回答 Srivatsan 在评论中的问题

_rvs是公共方法调用的特定于发行版的方法rvs. rvs是一个通用方法，它执行一些参数检查、添加位置和比例并设置属性self._size这是请求的随机变量数组的大小，然后调用分布特定方法._rvs或者它是通用的对应物。中的额外参数._rvs是形状参数，但由于本例中没有形状参数，*x and **y是多余的和未使用的。

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