Regresión sigmoidal con scipy, numpy, python, etc.

Tengo dos variables (x e y) que tienen una relación un tanto sigmoidal entre sí, y necesito encontrar algún tipo de ecuación de predicción que me permita predecir el valor de y, dado cualquier valor de x. Mi ecuación de predicción necesita mostrar la relación algo sigmoidal entre las dos variables. Por lo tanto, no puedo conformarme con una ecuación de regresión lineal que produce una línea. Necesito ver el cambio gradual y curvilíneo en la pendiente que se produce tanto a la derecha como a la izquierda de la gráfica de las dos variables.

Comencé a usar numpy.polyfit después de buscar en Google la regresión curvilínea y python, pero eso me dio los resultados terribles que puedes ver si ejecutas el código a continuación. ¿Alguien me puede mostrar cómo reescribir el código a continuación para obtener el tipo de ecuación de regresión sigmoidal que quiero?

Si ejecuta el código a continuación, puede ver que da una parábola orientada hacia abajo, que no es como debería ser la relación entre mis variables. En su lugar, debería haber más de una relación sigmoidal entre mis dos variables, pero con un ajuste perfecto con los datos que estoy usando en el código a continuación. Los datos en el código a continuación son los medios de un estudio de investigación de grandes muestras, por lo que tienen más poder estadístico del que podrían sugerir sus cinco puntos de datos. No tengo los datos reales del estudio de investigación de muestra grande, pero sí tengo los medios a continuación y sus desviaciones estándar (que no estoy mostrando). Preferiría simplemente trazar una función simple con los datos medios enumerados a continuación, pero el código podría volverse más complejo si la complejidad ofreciera mejoras sustanciales.

¿Cómo puedo cambiar mi código para mostrar el mejor ajuste de una función sigmoidal, preferiblemente usando scipy, numpy y python? Aquí está la versión actual de mi código, que debe ser corregida:

import numpy as np import matplotlib.pyplot as plt # Create numpy data arrays x = np.array([821,576,473,377,326]) y = np.array([255,235,208,166,157]) # Use polyfit and poly1d to create the regression equation z = np.polyfit(x, y, 3) p = np.poly1d(z) xp = np.linspace(100, 1600, 1500) pxp=p(xp) # Plot the results plt.plot(x, y, '.', xp, pxp, '-') plt.ylim(140,310) plt.xlabel('x') plt.ylabel('y') plt.grid(True) plt.show() 

EDITAR ABAJO: (Re-enmarcó la pregunta)

Su respuesta, y su velocidad, son muy impresionantes. Gracias, unutbu. Pero, para producir resultados más válidos, debo volver a enmarcar los valores de mis datos. Esto significa volver a emitir los valores de x como un porcentaje del valor máximo de x, mientras que volver a emitir los valores de y como un porcentaje de los valores de x en los datos originales. Intenté hacer esto con su código y se me ocurrió lo siguiente:

 import numpy as np import matplotlib.pyplot as plt import scipy.optimize # Create numpy data arrays ''' # Comment out original data #x = np.array([821,576,473,377,326]) #y = np.array([255,235,208,166,157]) ''' # Re-calculate x values as a percentage of the first (maximum) # original x value above x = np.array([1.000,0.702,0.576,0.459,0.397]) # Recalculate y values as a percentage of their respective x values # from original data above y = np.array([0.311,0.408,0.440,0.440,0.482]) def sigmoid(p,x): x0,y0,c,k=py = c / (1 + np.exp(-k*(x-x0))) + y0 return y def residuals(p,x,y): return y - sigmoid(p,x) p_guess=(600,200,100,0.01) (p, cov, infodict, mesg, ier)=scipy.optimize.leastsq(residuals,p_guess,args=(x,y),full_output=1,warning=True) ''' # comment out original xp to allow for better scaling of # new values #xp = np.linspace(100, 1600, 1500) ''' xp = np.linspace(0, 1.1, 1100) pxp=sigmoid(p,xp) x0,y0,c,k=p print('''\ x0 = {x0} y0 = {y0} c = {c} k = {k} '''.format(x0=x0,y0=y0,c=c,k=k)) # Plot the results plt.plot(x, y, '.', xp, pxp, '-') plt.ylim(0,1) plt.xlabel('x') plt.ylabel('y') plt.grid(True) plt.show() 

¿Me puede mostrar cómo arreglar este código revisado?
NOTA: Al volver a emitir los datos, esencialmente he girado el sigmoide 2d (x, y) sobre el eje z en 180 grados. Además, el 1.000 no es realmente un máximo de los valores de x. En cambio, 1.000 es una media del rango de valores de diferentes participantes de prueba en una condición de prueba máxima.


SEGUNDA EDICIÓN A CONTINUACIÓN:

Gracias, ubuntu. Leí cuidadosamente su código y busqué aspectos de él en la documentación de scipy. Ya que su nombre parece aparecer como un escritor de la documentación de scipy, espero que pueda responder las siguientes preguntas:

1.) ¿Lesssq () llama a residuales (), que luego devuelve la diferencia entre el vector y de entrada y el vector y devuelto por la función sigmoide ()? Si es así, ¿cómo explica la diferencia en las longitudes del vector y de entrada y del vector y devuelto por la función sigmoide ()?

2.) Parece que puedo llamar a leastsq () para cualquier ecuación matemática, siempre y cuando acceda a esa ecuación matemática a través de una función de residuos, que a su vez llama a la función matemática. ¿Es esto cierto?

3.) También, me doy cuenta de que p_guess tiene el mismo número de elementos que p. ¿Significa esto que los cuatro elementos de p_guess se corresponden en orden, respectivamente, con los valores devueltos por x0, y0, c y k?

4.) ¿La p que se envía como un argumento a las funciones residual () y sigmoid () es la misma p que se emitirá con leastsq (), y la función leastsq () usa esa p internamente antes de devolverla?

5.) ¿P y p_guess pueden tener cualquier número de elementos, dependiendo de la complejidad de la ecuación que se usa como modelo, siempre que el número de elementos en p sea igual al número de elementos en p_guess?

Utilizando scipy.optimize.leastsq :

 import numpy as np import matplotlib.pyplot as plt import scipy.optimize def sigmoid(p,x): x0,y0,c,k=p y = c / (1 + np.exp(-k*(x-x0))) + y0 return y def residuals(p,x,y): return y - sigmoid(p,x) def resize(arr,lower=0.0,upper=1.0): arr=arr.copy() if lower>upper: lower,upper=upper,lower arr -= arr.min() arr *= (upper-lower)/arr.max() arr += lower return arr # raw data x = np.array([821,576,473,377,326],dtype='float') y = np.array([255,235,208,166,157],dtype='float') x=resize(-x,lower=0.3) y=resize(y,lower=0.3) print(x) print(y) p_guess=(np.median(x),np.median(y),1.0,1.0) p, cov, infodict, mesg, ier = scipy.optimize.leastsq( residuals,p_guess,args=(x,y),full_output=1,warning=True) x0,y0,c,k=p print('''\ x0 = {x0} y0 = {y0} c = {c} k = {k} '''.format(x0=x0,y0=y0,c=c,k=k)) xp = np.linspace(0, 1.1, 1500) pxp=sigmoid(p,xp) # Plot the results plt.plot(x, y, '.', xp, pxp, '-') plt.xlabel('x') plt.ylabel('y',rotation='horizontal') plt.grid(True) plt.show() 

rendimientos

texto alternativo

con parametros sigmoides

 x0 = 0.826964424481 y0 = 0.151506745435 c = 0.848564826467 k = -9.54442292022 

Tenga en cuenta que para las versiones más nuevas de scipy (por ejemplo, 0.9) también existe la función scipy.optimize.curve_fit , que es más fácil de usar que la leastsq . Una discusión relevante sobre el ajuste de sigmoides usando curve_fit se puede encontrar aquí .

Edición: Se agregó una función de cambio de resize para que los datos sin procesar se puedan reajustar y cambiar para que se ajusten a cualquier cuadro delimitador deseado.


“Tu nombre parece aparecer como un escritor de la documentación de scipy”

DESCARGO DE RESPONSABILIDAD: No soy un escritor de documentación scipy. Solo soy un usuario, y un principiante en eso. Mucho de lo que sé acerca de leastsq viene de leer este tutorial , escrito por Travis Oliphant.

1.) ¿Lesssq () llama a residuales (), que luego devuelve la diferencia entre el vector y de entrada y el vector y devuelto por la función sigmoide ()?

¡Sí! exactamente.

Si es así, ¿cómo explica la diferencia en las longitudes del vector y de entrada y del vector y devuelto por la función sigmoide ()?

Las longitudes son las mismas:

 In [138]: x Out[138]: array([821, 576, 473, 377, 326]) In [139]: y Out[139]: array([255, 235, 208, 166, 157]) In [140]: p=(600,200,100,0.01) In [141]: sigmoid(p,x) Out[141]: array([ 290.11439268, 244.02863507, 221.92572521, 209.7088641 , 206.06539033]) 

Una de las cosas maravillosas de Numpy es que le permite escribir ecuaciones “vectoriales” que operan en arreglos completos.

 y = c / (1 + np.exp(-k*(x-x0))) + y0 

puede parecer que funciona en flotadores (de hecho lo haría) pero si haces x una matriz numpy, y c , k , x0 , y0 flota, entonces la ecuación define y es una matriz numpy de la misma forma que x . Entonces sigmoid(p,x) devuelve una matriz numpy. Hay una explicación más completa de cómo funciona esto en el numpybook (lectura requerida para usuarios serios de numpy).

2.) Parece que puedo llamar a leastsq () para cualquier ecuación matemática, siempre y cuando acceda a esa ecuación matemática a través de una función de residuos, que a su vez llama a la función matemática. ¿Es esto cierto?

Cierto. leastsq intenta minimizar la sum de los cuadrados de los residuos (diferencias). Busca en el espacio de parámetros (todos los valores posibles de p ) buscando la p que minimiza esa sum de cuadrados. Las x y y enviadas a los residuals son sus valores de datos sin procesar. Ellos son fijos. No cambian Es la p s (los parámetros en la función sigmoidea) lo que leastsq intenta minimizar.

3.) También, me doy cuenta de que p_guess tiene el mismo número de elementos que p. ¿Significa esto que los cuatro elementos de p_guess se corresponden en orden, respectivamente, con los valores devueltos por x0, y0, c y k?

¡Exacto así! Al igual que el método de Newton, leastsq necesita una conjetura inicial para p . Lo suministras como p_guess . Cuando veas

 scipy.optimize.leastsq(residuals,p_guess,args=(x,y)) 

se puede pensar que, como parte del algoritmo lesssq (en realidad, el algoritmo de Levenburg-Marquardt) como primer paso, lesssq llama a los residuals(p_guess,x,y) . Note la similitud visual entre

 (residuals,p_guess,args=(x,y)) 

y

 residuals(p_guess,x,y) 

Puede ayudarlo a recordar el orden y el significado de los argumentos en leastsq .

residuals , como sigmoid devuelve una matriz numpy. Los valores de la matriz se ajustan al cuadrado y luego se sumn. Este es el número a batir. p_guess se varía al menos. leastsq busca un conjunto de valores que minimice los residuals(p_guess,x,y) .

4.) ¿La p que se envía como un argumento a las funciones residual () y sigmoid () es la misma p que se emitirá con leastsq (), y la función leastsq () usa esa p internamente antes de devolverla?

Bueno no exactamente. Como ya se sabe, p_guess varía al menos leastsq busca el valor p que minimiza los residuals(p,x,y) . La p (er, p_guess ) que se envía a leastsq tiene la misma forma que la p que se devuelve con leastsq . Obviamente, los valores deberían ser diferentes a menos que seas un gran adivino 🙂

5.) ¿P y p_guess pueden tener cualquier número de elementos, dependiendo de la complejidad de la ecuación que se usa como modelo, siempre que el número de elementos en p sea igual al número de elementos en p_guess?

Sí. No he realizado pruebas de estrés con leastsq para una gran cantidad de parámetros, pero es una herramienta increíblemente poderosa.

No creo que vaya a obtener buenos resultados con un ajuste polinomial de ningún grado, ya que todos los polinomios van al infinito para una X suficientemente grande y pequeña, pero una curva sigmoidea se aproximará asintóticamente a algún valor finito en cada dirección.

No soy un progtwigdor de Python, así que no sé si numpy tiene una rutina de ajuste de curvas más general. Si tiene que hacer su propio rollo, quizás este artículo sobre la regresión logística le dará algunas ideas.

Para la regresión logística en Python, scikits-learn expone el código de adaptación de alto rendimiento:

http://scikit-learn.sourceforge.net/modules/linear_model.html#logistic-regression

Como lo señaló @unutbu anterior, scipy ahora ofrece scipy.optimize.curve_fit que poseen una llamada menos complicada. Si alguien quiere una versión rápida de cómo se vería el mismo proceso en esos términos, presento un ejemplo mínimo a continuación:

 def sigmoid(x, k, x0): return 1.0 / (1 + np.exp(-k * (x - x0))) # Parameters of the true function n_samples = 1000 true_x0 = 15 true_k = 1.5 sigma = 0.2 # Build the true function and add some noise x = np.linspace(0, 30, num=n_samples) y = sigmoid(x, k=true_k, x0=true_x0) y_with_noise = y + sigma * np.random.randn(n_samples) # Sample the data from the real function (this will be your data) some_points = np.random.choice(1000, size=30) # take 30 data points xdata = x[some_points] ydata = y_with_noise[some_points] # Fit the curve popt, pcov = curve_fit(return_sigmoid, xdata, ydata) estimated_k, estimated_x0 = popt # Plot the fitted curve y_fitted = sigmoid(x, k=estimated_k, x0=estimated_x0) # Plot everything for illustration fig = plt.figure() ax = fig.add_subplot(111) ax.plot(x, y_fitted, '--', label='fitted') ax.plot(x, y, '-', label='true') ax.plot(xdata, ydata, 'o', label='samples') ax.legend() 

El resultado de esto se muestra en la siguiente figura:

introduzca la descripción de la imagen aquí