library - scipy python 3
Regresión sigmoidal con scipy, numpy, python, etc. (4)
Tengo dos variables (x e y) que tienen una relación algo 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. ¿Puede alguien mostrarme cómo reescribir el código a continuación para obtener el tipo de ecuación de regresión sigmoidal que deseo?
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, necesito 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=p
y = 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 residuals (), 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?
Como señaló @unutbu anterior, scipy ahora ofrece scipy.optimize.curve_fit que posee 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:
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 programador 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
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
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 utilizando curve_fit se puede encontrar here .
Edición: Se agregó una función de cambio de resize para que los datos sin procesar se puedan reescala y se desplacen para que se ajusten a cualquier cuadro de límite 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 residuals (), 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 para ser 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 suma 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 suma de cuadrados. Las x y y enviadas a los residuals son sus valores de datos sin procesar. Ellos son fijos. No cambian Son las 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 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 suman. Este es el número a vencer. 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 devuelve 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 probado los leastsq de leastsq para una gran cantidad de parámetros, pero es una herramienta increíblemente poderosa.