在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方法的形狀有多好。這些分佈是爲單變量分佈而設計的，可能不適用於多變量分佈情況，或者可能需要一些重構。