Uso de numpy.optimize.leastsq() o cómo ahorrar trabajo.

Como sospechaba en la entrada de ayer, es posible ‘engañar’ a la función optimize.leastsq de numpy para trabajar con arrays n-dimensionales.
El resultado es mejor (por poco😉 ) al obtenido mediante el algoritmo de Gauss-Newton ayer, aunque no me cabe la menor duda de que será mucho más eficiente.

A trabajar.

Para que esta función haga su trabajo necesitaremos darle algunas cosas hechas:
  • Una función que calcule la imagen (matemática) de nuestra función. La llamaré
  • Una función que calcule el residual. La llamaré residual.
  • Un arras (no es obligatorio pero sí más cómodo) que contenga las variables a optimizar. Lo llamaré beta.
Con esto podremos poner a trabajar a nuestra función de la siguiente forma:

def GaussianFunct( beta, x, y):
 H,B,sigma,x0,y0 = beta

 exponente = -( (x-x0)**2 +(y-y0)**2) / (2*sigma*sigma)
 return H*np.e**(exponente) + B

def residual(beta, x,y,z):
 return (z-GaussianFunct(beta,x,y))**2

Ahora vamos con la optimización:

beta,cov,infodict,mesg,ier = optimize.leastsq( residual, beta, args=(meshX, meshY, Z), full_output=True)

Al tratarse mi caso de una función tridimensional (x,y,z) (no en el ámbito matemático, pero así nos entendemos) es cómodo evaluarla mediante ‘mallas’ obtenidas mediante meshgrid(), pero optimize.leastsq() no los acepta. No hay problema, las aplanamos primero mediante ravel() o reshape():

z = star[0:24,0:24]
x = np.arange(24)
y = np.arange(24)

meshx,meshy = np.meshgrid(x,y)

Z = np.ravel( z )
meshX = np.ravel( meshx )
meshY = np.ravel( meshy )

beta = np.array([ 100 ,6, 1.6,12, 12])

Resultados

El resultado es muy similar al de la entrada anterior en términos de error máximo porcentual, pero mirando con un poco más de detenimiento la distribución del error en la superficie no es la misma…

errorNumpy

Dejo aquí el notebook para que se vea más claro el flujo del programa😉.

Esta entrada fue publicada en matemáticas, python y etiquetada , , , , . Guarda el enlace permanente.

10 respuestas a Uso de numpy.optimize.leastsq() o cómo ahorrar trabajo.

  1. Juanlu001 dijo:

    Te voy a fichar de colaborador en el blog…😛

  2. Juanlu001 dijo:

    De hecho, creo que aún podrías ahorrar más trabajo usando el módulo `scipy.stats`:

    http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.rv_continuous.fit.html

    Las distribuciones tienen su propio método `fit`, que utiliza el método de máxima verosimilitud y resuelve el problema de optimización usando el algoritmo “downhill simplex” (aunque esto se puede cambiar).

    Pensaba tratar el tema en el blog (como continuación a esto (http://pybonacci.wordpress.com/2012/04/21/estadistica-en-python-con-scipy/) pero hombre, si te me adelantas podríamos llegar a un acuerdo jeje😛

    • Samuel dijo:

      ¿Sugieres que mediante una distribución estadística podríamos ajustar mejor los datos? ¿qué ventaja puede tener eso si conozco a priori la función que mejor resultado va a dar? ¿sólo la velocidad?

      • Juanlu001 dijo:

        O sea, en realidad ya estás usando una distribución estadística. Este sería otro método distinto: tú ya estás optimizando bastante, así que no sé realmente si ganarías mucho o poco. Pero vamos, que en teoría se puede hacer de la otra manera también.

  3. Samuel dijo:

    He dejado al final de la entrada un enlace al notebook de python impreso, para quien tenga dudas!

  4. Juanlu001 dijo:

    Hm, para el notebook lo mejor es que subas el archivo .ipynb a algún sitio (File -> Download as… -> .ipynb) y que después le metas esa URL aquí:

    http://nbviewer.ipython.org/

    por ejemplo, aquí tengo subido un notebook:

    https://raw.github.com/Pybonacci/poliastro/master/examples/Going%20to%20Mars%20with%20Python%20using%20poliastro.ipynb

    Y si meto la URL ahí:

    http://nbviewer.ipython.org/urls/raw.github.com/Pybonacci/poliastro/master/examples/Going%2520to%2520Mars%2520with%2520Python%2520using%2520poliastro.ipynb

    Ya lo tengo🙂

Responder

Introduce tus datos o haz clic en un icono para iniciar sesión:

Logo de WordPress.com

Estás comentando usando tu cuenta de WordPress.com. Cerrar sesión / Cambiar )

Imagen de Twitter

Estás comentando usando tu cuenta de Twitter. Cerrar sesión / Cambiar )

Foto de Facebook

Estás comentando usando tu cuenta de Facebook. Cerrar sesión / Cambiar )

Google+ photo

Estás comentando usando tu cuenta de Google+. Cerrar sesión / Cambiar )

Conectando a %s