在 scipy 中创建新的发行版

Posted 2023-03-15

【中文标题】在 scipy 中创建新的发行版

【英文标题】：Creating new distributions in scipy

【发布时间】：2012-05-27 13:23:12

【问题描述】：

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

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（最近一次调用最后一次）：文件“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 行，在 调用 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-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) and f(b) must have different sign

编辑

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

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)

【问题讨论】：

那么当我们执行上述操作并调用randoms = getDistribution(Mydata) 然后调用randoms = randoms.rvs(size=1000) 时，它会在类中执行三个def 吗？即计算pdf，整合它等？

我确实让我的随机数遵循数据分布，但我想对其进行平滑处理，使其不完全遵循数据分布。我一直在手动调整kernel 的带宽来做到这一点。例如，我们如何指定 PDF 函数，然后使用 PDF 函数使用 Metropolis Hastings 创建随机数。

【参考方案1】：

特别是您的回溯：

rvs 使用 cdf 的倒数 ppf 创建随机数。由于您没有指定 ppf，它是通过寻根算法brentq 计算的。 brentq 使用下限和上限来搜索函数为零的值 at （查找 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 有一个 resample 方法。可以通过从数据中采样并添加高斯噪声来生成随机样本。因此，这将比使用 ppf 方法的通用 rvs 更快。我会写一个 ._rvs 方法，它只调用 gaussian_kde 的 resample 方法。

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

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

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

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

【讨论】：

太棒了，谢谢，这很有帮助。是否有某种方法可以使用与 scipy 相同的方法从 cdf 预先计算 ppf，从而提高效率？我注意到每次 rv() 调用都会调用 _cdf() 很多次。

我在 rvs 和 ppf 上添加了一些 cmets。还有一条评论：如果您有带有肥尾的数据，那么 gaussian_kde 的尾巴就不会很好。当我考虑编写一个类似的分布子类时，我会尝试使用帕累托尾巴。我在一个论坛上读到了关于这个的评论，matlab 有一个帕累托尾分布。

@user333700 我想知道问题中的编辑发布的答案中的 def _rvs(self, *x, **y) 有什么作用？？