Creazione di nuove distribuzioni in scipy

Question

Sto provando a creare una distribuzione basata su alcuni dati che ho, quindi disegno a caso da quella distribuzione. Ecco quello che ho:

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 

Penso che questo stia facendo quello che voglio, ma ho spesso ottenere un errore (vedi sotto), quando cerco di fare d.rvs(), e d.rvs(100) non funziona mai. Sto facendo qualcosa di sbagliato? C'è un modo più semplice o migliore per farlo? Se si tratta di un bug in scipy, c'è un modo per aggirarlo?

Infine, c'è più documentazione sulla creazione di distribuzioni personalizzate da qualche parte? Il meglio che ho trovato è la documentazione scipy.stats.rv_continuous, che è piuttosto spartana e non contiene esempi utili.

Il 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

Modifica

Per chi fosse curioso, seguendo i consigli nella risposta di seguito, ecco il codice che funziona:

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

particolare al vostro traceback:

rvs utilizza la i nverse del cdf, ppf, per creare numeri casuali. Poiché non stai specificando ppf, viene calcolato da un algoritmo di rootfinding, brentq. brentq utilizza i limiti inferiore e superiore su dove deve cercare il valore at con la funzione è zero (trova x tale che cdf (x) = q, q è quantile).

L'impostazione predefinita per i limiti, xa e xb, è troppo piccola nell'esempio. Le seguenti opere per me con SciPy 0.9.0, xa, xb possono essere impostati quando si crea l'istanza funzione di

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) 

Attualmente v'è una richiesta di pull per SciPy per migliorare questo, quindi nella prossima release xa e xb sarà essere espanso automaticamente per evitare l'eccezione f(a) and f(b) must have different signs.

Non c'è molta documentazione su questo, il più semplice è seguire alcuni esempi (e chiedere sulla mailing list).

edit: Oltre

pdf: Dal momento che hai la funzione di densità data anche da gaussian_kde, vorrei aggiungere il metodo _pdf, che renderà alcuni calcoli più efficiente.

EDIT2: Oltre

rvs: Se siete interessati a generare numeri casuali, quindi gaussian_kde ha un metodo di ricampionamento. Campioni casuali possono essere generati campionando dai dati e aggiungendo rumore gaussiano. Quindi, questo sarà più veloce rispetto ai generici rvs che usano il metodo ppf. Scriverò un metodo ._rvs che chiama semplicemente il metodo resample di gaussian_kde.

precomputing ppf: Non conosco alcun modo generale per precomputare il ppf.Tuttavia, il modo in cui ho pensato di farlo (ma mai provato finora) è precomputare il ppf in molti punti e quindi utilizzare l'interpolazione lineare per approssimare la funzione ppf.

Edit3: circa _rvs per rispondere alla domanda di Srivatsan nel commento

_rvs è il metodo specifico di distribuzione che viene chiamato dal metodo pubblico rvs. rvs è un metodo generico che esegue il controllo degli argomenti, aggiunge la posizione e la scala e imposta l'attributo self._size, che corrisponde alla dimensione dell'array richiesto di variabili casuali, quindi chiama il metodo specifico di distribuzione ._rvs o la sua controparte generica. Gli argomenti aggiuntivi in ​​._rvs sono parametri di forma, ma poiché non vi sono in questo caso, *x e **y sono ridondanti e non utilizzati.

Non so quanto bene il size o la forma del metodo .rvs funzioni nel caso multivariato. Queste distribuzioni sono progettate per distribuzioni univariate e potrebbero non funzionare completamente per il caso multivariato o potrebbero richiedere alcuni rimodellamenti.