Erstellen neuer Distributionen in scipy

Question

Ich versuche, eine Verteilung basierend auf einigen Daten zu erstellen, die ich habe, und dann zufällig aus dieser Verteilung zu zeichnen. Hier ist, was ich habe:

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 

Ich denke, das tut, was ich will, aber ich bekommen häufig einen Fehler (siehe unten), wenn ich versuche d.rvs() zu tun, und d.rvs(100) nie funktioniert. Mache ich etwas falsch? Gibt es einen leichteren oder besseren Weg dies zu tun? Wenn es ein Fehler in scipy ist, gibt es einen Weg, um es zu umgehen?

Schließlich gibt es weitere Dokumentation zum Erstellen von benutzerdefinierten Distributionen irgendwo? Das Beste, was ich gefunden habe, ist die Dokumentation scipy.stats.rv_continuous, die ziemlich spartanisch ist und keine nützlichen Beispiele enthält.

Die Zurückverfolgungs:

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

bearbeiten

Für diejenigen neugierig, im Anschluss an die Beratung in der Antwort unten, hier Code ist, das funktioniert:

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

speziell auf Ihre Traceback:

rvs verwendet die i nverse des cdf, ppf, um Zufallszahlen zu erstellen. Da Sie ppf nicht angeben, wird es von einem Rootfinding-Algorithmus brentq berechnet. brentq verwendet untere und obere Grenzen auf, wo es nach dem Wert bei mit der Funktion suchen sollte, ist Null (find x so, dass cdf (x) = q, q ist Quantil).

Die Standardwerte für die Grenzwerte xa und xb sind in Ihrem Beispiel zu klein. Die folgenden Werke für mich mit scipy 0.9.0, xa können xb eingestellt werden, wenn die Funktion Instanz

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) 

Es gibt derzeit eine Pull-Anforderung für scipy dies zu verbessern, so dass in der nächsten Version xa und xb wird automatisch erweitert werden, um die Ausnahme f(a) and f(b) must have different signs zu vermeiden.

Es gibt nicht viel Dokumentation darüber, am einfachsten ist es, ein paar Beispiele zu folgen (und auf der Mailingliste zu fragen).

bearbeiten: Neben

pdf: Da Sie die Dichtefunktion auch durch gaussian_kde gegeben, würde ich die _pdf Methode, fügen die einige Berechnungen effizienter machen.

edit2: Zusatz

rvs: Wenn Sie bei der Erzeugung von Zufallszahlen interessiert sind, dann hat gaussian_kde eine Resampling-Methode. Random Samples können durch Abtasten der Daten und Hinzufügen von Gauß-Rauschen erzeugt werden. Also, dies wird schneller sein als die generischen rvs mit der ppf-Methode. Ich würde eine ._rvs-Methode schreiben, die nur Gaussian_kde's Resample-Methode aufruft.

vorberechnen ppf: Ich kenne keine allgemeine Möglichkeit, die ppf vorberechnen.Die Art und Weise, wie ich daran gedacht habe (aber bisher noch nie versucht habe) ist, den ppf an vielen Punkten vorzuberechnen und dann lineare Interpolation zu verwenden, um die ppf-Funktion zu approximieren.

EDIT3: über _rvs zu beantworten Srivatsan Frage im Kommentar

_rvs ist die Verteilung spezifische Methode, die rvs von der öffentlichen Methode aufgerufen wird. rvs ist eine generische Methode, die einige Argumente überprüft, fügt Ort und Skalierung hinzu und legt das Attribut self._size fest, das die Größe des angeforderten Arrays zufälliger Variablen darstellt, und ruft dann die verteilungsspezifische Methode ._rvs oder das generische Gegenstück auf. Die zusätzlichen Argumente in ._rvs sind Formparameter, aber da es in diesem Fall keine gibt, sind *x und **y redundant und unbenutzt.

Ich weiß nicht, wie gut die size oder die Form der .rvs Methode im multivariaten Fall funktioniert. Diese Verteilungen sind für univariate Verteilungen konzipiert und funktionieren möglicherweise nicht vollständig für den multivariaten Fall oder benötigen möglicherweise einige Umformungen.