首页 文章

在scipy中创建新的发行版

提问于
浏览 次
11

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

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（最近一次调用最后一次）：文件“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-p ackages / 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）必须有不同的符号

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)
python numpy scipy

1 回答

  • 7

    特别是你的追溯：

    rvs使用cdf的倒数ppf来创建随机数 . 由于您没有指定ppf，因此它由rootfinding算法 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请求来改进这一点，因此在下一个版本中 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方法 .

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

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

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

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

    回复于

相关问题