Creando nuevas distribuciones en scipy

Question

Estoy tratando de crear una distribución basada en algunos datos que tengo, luego dibujo al azar de esa distribución. Aquí es lo que tengo:

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 

creo que esto es hacer lo que yo quiero que, con frecuencia, pero me sale un error (véase más adelante) cuando intento hacer d.rvs(), y d.rvs(100) nunca funciona. ¿Estoy haciendo algo mal? ¿Hay una manera más fácil o mejor de hacer esto? Si es un error en scipy, ¿hay alguna forma de evitarlo?

Finalmente, ¿hay más documentación sobre cómo crear distribuciones personalizadas en alguna parte? Lo mejor que he encontrado es la documentación scipy.stats.rv_continuous, que es bastante espartana y no contiene ejemplos útiles.

El rastreo:

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

Editar

Para los curiosos, siguiendo el consejo de la respuesta a continuación, aquí está el código que funciona:

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

específicamente a su rastreo:

rvs usa el i nverse del cdf, ppf, para crear números aleatorios. Como no especifica ppf, se calcula mediante un algoritmo de búsqueda de root, brentq. brentq usa los límites inferior y superior en donde debe buscar el valor en con la función es cero (encuentre x tal que cdf (x) = q, q es cuantil).

El valor predeterminado para los límites, xa y xb, son demasiado pequeños en su ejemplo. Las siguientes obras para mí con scipy 0.9.0, xa, xb se pueden establecer cuando se crea la instancia de función

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) 

Actualmente hay una solicitud de extracción de scipy para mejorar esto, así que en la próxima versión xa y se xb se expandirá automáticamente para evitar la excepción f(a) and f(b) must have different signs.

No hay mucha documentación sobre esto, lo más fácil es seguir algunos ejemplos (y preguntar en la lista de correo).

edición: Además

pdf: Puesto que usted tiene la función de densidad dada también por gaussian_kde, yo añadiría el método _pdf, lo que hará que algunos cálculos más eficiente.

Edit2: Además

rvs: Si usted está interesado en la generación de números aleatorios, entonces gaussian_kde tiene un método de remuestreo. Las muestras aleatorias pueden generarse al tomar muestras de los datos y agregar ruido gaussiano. Entonces, esto será más rápido que los rvs genéricos usando el método ppf. Escribiría un método ._rvs que simplemente llama al método de remuestreo de gaussian_kde.

precomputing ppf: No conozco ninguna forma general de predecir el ppf.Sin embargo, la forma en que pensé hacerlo (pero nunca lo intenté) es calcular previamente la ppf en muchos puntos y luego usar la interpolación lineal para aproximar la función ppf.

Edit3: sobre _rvs para responder a la pregunta de Srivatsan en el comentario

_rvs es el método específico de distribución citado por el método público rvs. rvs es un método genérico que verifica algunos argumentos, agrega ubicación y escala, y establece el atributo self._size que es el tamaño de la matriz solicitada de variables aleatorias, y luego llama al método específico de distribución ._rvs o su contraparte genérica. Los argumentos adicionales en ._rvs son parámetros de forma, pero como no los hay, en este caso, *x y **y son redundantes y no se utilizan.

No sé qué tan bien funciona el size o la forma del método .rvs en el caso de multivariante. Estas distribuciones están diseñadas para distribuciones univariadas, y es posible que no funcionen completamente para el caso multivariable, o que necesiten algunas remodelaciones.