Inicio Desarrollo Regresión no lineal en Pyton

Regresión no lineal en Pyton

45684

Importación de los paquetes básicos

[sourcecode language=»python» wraplines=»false» collapse=»false»]
import numpy as np # paquete básico para calculos científicos
import matplotlib.pyplot as plt # para crear graficos
%matplotlib inline # muestra los graficos en la linea de comandos de Python
[/sourcecode]

Creación de diferentes conjuntos de datos

Función cúbica ( $y = a x^{3} + b x^{2} + c x + d$ )

[sourcecode language=»python» wraplines=»false» collapse=»false» firstline=»4″]
x = np.arange(-5.0, 5.0, 0.1)
##Puede ajustar la pendiente y la intersección para verificar los cambios del gráfico
y = 1*(x**3) + 1*(x**2) + 1*x + 3
y_noise = 20 * np.random.normal(size=x.size)
ydata = y + y_noise
plt.plot(x, ydata, ‘bo’)
plt.plot(x,y, ‘r’)
plt.ylabel(‘Variable dependiente’)
plt.xlabel(‘Variable indepdendiente’)
plt.show()
[/sourcecode]

Función cuadrática ( $y = x^{2}$ )

[sourcecode language=»python» wraplines=»false» collapse=»false» firstline=»14″]
x = np.arange(-5.0, 5.0, 0.1)
##Se puede ajustar la pendiente y la intersección para verificar los cambios en el gráfico
y = np.power(x,2)
y_noise = 2 * np.random.normal(size=x.size)
ydata = y + y_noise
plt.plot(x, ydata, ‘bo’)
plt.plot(x,y, ‘r’)
plt.ylabel(‘Variable dependiente’)
plt.xlabel(‘Variable indepdiendente’)
plt.show()
[/sourcecode]

Exponencial ( $y = a + b c^{x}$ )

[sourcecode language=»python» wraplines=»false» collapse=»false» firstline=»24″]
X = np.arange(-5.0, 5.0, 0.1)
##Se puede ajustar la pendiente y la intersección para verificar los cambios en el gráfico
Y= np.exp(X)
plt.plot(X,Y)
plt.ylabel(‘Variable Dependiente’)
plt.xlabel(‘Variable Independiente’)
plt.show()
[/sourcecode]

Logarítmica ( $y = \log(x)$ )

[sourcecode language=»python» wraplines=»false» collapse=»false» firstline=»31″]
X = np.arange(-5.0, 5.0, 0.1)
Y = np.log(X)
plt.plot(X,Y)
plt.ylabel(‘Variable Dependiente’)
plt.xlabel(‘Variable Independiente’)
plt.show()
[/sourcecode]

Sigmoidal/Logística ( $y = a + \dfrac{b}{1 + c^{(x - d)}}$ )

[sourcecode language=»python» wraplines=»false» collapse=»false» firstline=»37″]
X = np.arange(-5.0, 5.0, 0.1)
Y = 1-4/(1+np.power(3, X-2))
plt.plot(X,Y)
plt.ylabel(‘Variable Dependiente’)
plt.xlabel(‘Variable Independiente’)
plt.show()
[/sourcecode]

Ejemplo completo

Descarga del fichero de datos

[sourcecode language=»python» wraplines=»false» collapse=»false» firstline=»1″]
import numpy as np
import pandas as pd
#downloading dataset
!wget -nv -O china_gdp.csv https://s3-api.us-geo.objectstorage.softlayer.net/cf-courses-data/CognitiveClass/ML0101ENv3/labs/china_gdp.csv
df = pd.read_csv("china_gdp.csv")
df.head(10) # mostramos los primeros datos
[/sourcecode]

	Year	Value
0	1960	5.918412e+10
1	1961	4.955705e+10
2	1962	4.668518e+10
3	1963	5.009730e+10
4	1964	5.906225e+10
5	1965	6.970915e+10
6	1966	7.587943e+10
7	1967	7.205703e+10
8	1968	6.999350e+10
9	1969	7.871882e+10

Mostrado de los datos

[sourcecode language=»python» wraplines=»false» collapse=»false» firstline=»7″]
plt.figure(figsize=(8,5))
x_data, y_data = (df["Year"].values, df["Value"].values)
plt.plot(x_data, y_data, ‘ro’)
plt.ylabel(‘GDP’)
plt.xlabel(‘Year’)
plt.show()
[/sourcecode]

Eligiendo un modelo

A primera vista, determinamos que la función lógica podría ser una buena primera aproximación, ya que tiene la propiedad de comenzar con un crecimiento leve, aumentando en el medio y luego descendiendo nuevamente hacia el final; como vimos anteriormente y simplificando, la formula es $y = \dfrac{1}{1 + e^{\beta_1(x - \beta_2)}}$ donde $\beta_1$ controla lo llano de la curva y $\beta_2$ lleva la curva sobre el eje x.

[sourcecode language=»python» wraplines=»false» collapse=»false» firstline=»13″]
#define el modelo
def sigmoid(x, Beta_1, Beta_2):
y = 1 / (1 + np.exp(-Beta_1*(x-Beta_2)))
return y

#ejecutamos un primer modelo
beta_1 = 0.10
beta_2 = 1990.0

#función logística
Y_pred = sigmoid(x_data, beta_1 , beta_2)

#predicción de puntos
plt.plot(x_data, Y_pred*15000000000000.)
plt.plot(x_data, y_data, ‘ro’)
[/sourcecode]

Nuestra tarea aquí es encontrar los mejores parámetros para nuestro modelo. Normalicemos primero nuestro $x$ e $y$ :

[sourcecode language=»python» wraplines=»false» collapse=»false» firstline=»28″]
xdata =x_data/max(x_data)
ydata =y_data/max(y_data)
[/sourcecode]

Para buscar los mejores parámetros podemos utilizar curve_fit la cual utiliza cuadrados mínimos no lineales para cuadrar con la función sigmoide. popt son nuestros parámetros optimizados.

[sourcecode language=»python» wraplines=»false» collapse=»false» firstline=»30″]
from scipy.optimize import curve_fit
popt, pcov = curve_fit(sigmoid, xdata, ydata)
#imprimir los parámetros finales
print(" beta_1 = %f, beta_2 = %f" % (popt[0], popt[1]))
[/sourcecode]

beta_1 = 690.447530, beta_2 = 0.997207

Dibujamos el modelo

[sourcecode language=»python» wraplines=»false» collapse=»false» firstline=»34″]
x = np.linspace(1960, 2015, 55)
x = x/max(x)
plt.figure(figsize=(8,5))
y = sigmoid(x, *popt)
plt.plot(xdata, ydata, ‘ro’, label=’data’)
plt.plot(x,y, linewidth=3.0, label=’fit’)
plt.legend(loc=’best’)
plt.ylabel(‘GDP’)
plt.xlabel(‘Year’)
plt.show()
[/sourcecode]

Vemos la exactitud de nuestro modelo

[sourcecode language=»python» wraplines=»false» collapse=»false» firstline=»44″]
# divide los datos en entrenamiento y prueba
msk = np.random.rand(len(df)) < 0.8
train_x = xdata[msk]
test_x = xdata[~msk]
train_y = ydata[msk]
test_y = ydata[~msk]

# construye el modelo utilizando el set de entrenamiento
popt, pcov = curve_fit(sigmoid, train_x, train_y)

# predecir utilizando el set de prueba
y_hat = sigmoid(test_x, *popt)

# evaluation
print("Promedio de error absoluto: %.2f" % np.mean(np.absolute(y_hat – test_y)))
print("Suma residual de cuadrados (MSE): %.2f" % np.mean((y_hat – test_y) ** 2))
from sklearn.metrics import r2_score
print("R2-score: %.2f" % r2_score(y_hat , test_y) )
[/sourcecode]

Promedio de error absoluto: 0.04
Suma residual de cuadrados (MSE): 0.00
R2-score: 0.97