Tworzenie nowych dystrybucji w scipy

Question

Próbuję utworzyć dystrybucję na podstawie niektórych danych, które mam, a następnie losowo losować z tej dystrybucji. Oto co mam:

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 

myślę, że robi to, co chcę go, ale często pojawia się błąd (patrz niżej), gdy próbuję zrobić d.rvs() i d.rvs(100) nigdy nie działa. czy robię coś źle? Czy istnieje lepszy lub lepszy sposób na zrobienie tego? Jeśli jest to błąd w scipie, czy jest jakiś sposób na obejście tego?

Wreszcie, czy jest gdzieś więcej dokumentacji na temat tworzenia niestandardowych dystrybucji? Najlepsze, co znalazłem, to dokumentacja scipy.stats.rv_continuous, która jest dość spartańska i nie zawiera użytecznych przykładów.

traceback:

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

Edit

Dla tych, nowoczesny, po poradę w poniższej odpowiedzi, oto kod, który działa:

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

szczególności do traceback:

rvs używa i nverse z cdf, ppf, aby utworzyć losowe liczby. Ponieważ nie określasz ppf, jest on obliczany za pomocą algorytmu rootfinding, brentq. brentq używa dolnych i górnych granic, w których powinien szukać wartości przy funkcji jest zero (znaleźć x takie, że cdf (x) = q, q jest kwantylem).

Domyślne wartości limitów: xa i xb są w tym przykładzie zbyt małe. Następujące prace dla mnie z scipy 0.9.0, xa, xb można ustawić podczas tworzenia instancji funkcji

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) 

Obecnie wniosek ciągnąć za scipy poprawić to, więc w następnym wydaniu xa i xb będzie być automatycznie rozszerzany, aby uniknąć wyjątku od f(a) and f(b) must have different signs.

Nie ma wiele dokumentacji na ten temat, najłatwiej jest wykonać kilka przykładów (i zapytaj na listę mailingową).

edit: dodatek

pdf: Skoro masz funkcję gęstości również podane przez gaussian_kde, chciałbym dodać metodę _pdf, która sprawi, że niektóre obliczenia bardziej wydajne.

Edit2: dodatek

kupię: Jeśli jesteś zainteresowany generowanie liczb losowych, a następnie gaussian_kde ma metodę resample. Losowe próbki mogą być generowane przez próbkowanie z danych i dodawanie szumu gaussowskiego. Tak więc będzie to szybsze niż ogólne rvs przy użyciu metody ppf. Chciałbym napisać metodę ._rvs, która po prostu wywołuje metodę resample do gaussian_kde.

precomputing ppf: Nie znam żadnego ogólnego sposobu wstępnego obliczenia ppf.Jednak sposób, w jaki myślałem o zrobieniu tego (ale nigdy dotąd nie próbowałem) polega na wstępnym obliczeniu ppf w wielu punktach, a następnie zastosowaniu interpolacji liniowej w celu przybliżenia funkcji ppf.

Edit3: około _rvs odpowiedzieć na pytanie Srivatsan w komentarzu

_rvs jest specyficzna metoda dystrybucji, które nazywa się metodą publicznego rvs. rvs to ogólna metoda, która wykonuje pewne sprawdzanie argumentów, dodaje położenie i skalę oraz ustawia atrybut self._size, który jest rozmiarem żądanej tablicy zmiennych losowych, a następnie wywołuje metodę specyficzną dla dystrybucji ._rvs lub jej ogólny odpowiednik. Dodatkowe argumenty w parametrach ._rvs są parametrami kształtu, ale ponieważ nie ma w tym przypadku żadnych, *x i **y są nadmiarowe i nieużywane. Nie mam pojęcia, jak dobrze size lub kształt metody .rvs działa w przypadku wielowymiarowym. Dystrybucje te są przeznaczone dla dystrybucji jednowymiarowych i mogą nie działać w pełni dla przypadku wielowymiarowego lub mogą wymagać pewnych zmian.