En esta sección trabajaremos con un dataset conocido como House Sales in King County, USA, presentado en la plataforma de Kaggle [4], el cual es un gran dataset para evaluar simples modelos de regresión. Los registros contienen distintas características asociadas a las ventas de casas en la localidad King County, entre mayo de 2014 y mayo de 2015, las cuales vienen descritas en el dataset, como la cantidad de habitaciones, cantidad de baños, número de pisos, etc. Donde una de las variables a estudiar corresponde al precio en el cual se vendió la casa.
a) Construya un dataframe con los datos a analizar descargándolos desde la plataforma como se indicó. Explique por qué se realiza la línea 4.
In [1]:
import pandas as pd
import numpy as np
df = pd.read_csv("kc_house_data.csv")
df.drop(['id','date','zipcode',],axis=1,inplace=True)
df.head()
Out[1]:
La línea 4 se encarga de eliminar parámetros que no son importantes para la valoración de la vivienda. La función drop toma un arreglo de atributos y los elimina del dataset.
b) Describa brevemente el dataset a utilizar.
In [2]:
df.shape
df.info()
df.describe()
Out[2]:
Como se mustra en la salida anterior, el dataset corresponde a la información de 21613 casas, en donde cada una tiene 18 atributos asociados. Estos atributos corresponden a precio, número de baños, número de dormitorios, superficie del living, superficie total del terreno, número de pisos, entre otros más especificados en las columnas. Además, se presentan la media, la desviaciñon estándar y los cuartiles de cada atributo sobre el dataset total.
c) Normalice los datos antes de trabajar y aplique una transformación adecuada a la variable a predecir. Explique la importancia/conveniencia de realizar estas dos operaciones.
In [3]:
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
df_scaled = pd.DataFrame(scaler.fit_transform(df), columns=df.columns)
df_scaled['price'] = np.log(df['price'])
La normalización es necesaria debido a que la mayoría de los algoritmos clasificadores hacen uso de la distancia euclidiana para calcular la distancia entre dos puntos en el dataset. Cuando existe un atributo con un rango de valores mucho mayor en comparación a los demás, este gobernará el calculo de la distancia. Esto puede producir, por ejemplo, que el algoritmo de gradiente (SGD) converja mucho más lento en comparación a una instancia normalizada.
La transformación sobre el precio, por otro lado, es necesaria para utilizar regresión lineal. Las regresiones lineales, tal y como dice su nombre, modelan utilizando rectas. Estas rectas luego son comparadas con los valores del precio para probar que tan bien se ajustan a nuestros valores. No tiene sentido realizar esta comapracion si el precio no sigue una tendencia lineal. Por ejemplo, se estaría intentando ajustar una recta sobre una función cuadratica, lo cual no es correcto. Al aplicar el logaritmo, nos aseguramos que el precio, nuestro y , sea siempre lineal.
d) Realice una regresión lineal de mínimos cuadrados básica. Explique la importancia/conveniencia del paso 4 y los argumentos que se deben entregar a la función que implementa la regresión lineal.
In [4]:
import sklearn.linear_model as lm
X = df_scaled.iloc[:,1:] #use .ix instead, in older pandas version
N = X.shape[0]
X.insert(X.shape[1], 'intercept', np.ones(N))
y = df_scaled['price']
#mascara estatica con el 70% de los datos
mascara = np.zeros(len(X))
limit = int(len(X)*0.7)
mascara[:limit] = 1
istrain = mascara== 1
Xtrain = X[istrain]
ytrain = y[istrain]
Xtest = X[np.logical_not(istrain)]
ytest = y[np.logical_not(istrain)]
linreg = lm.LinearRegression(fit_intercept = False)
linreg.fit(Xtrain, ytrain)
Out[4]:
En el paso 4, se incorpora una columna de unos para representar a la constante $\beta_0$ en el modelo de la regresión lineal, que corresponde a nuestro intercepto. Esta constante proviene de la expresión:
$H = \beta^TX + \beta_0$
en donde $\beta$ es la matriz de los valores libres asociados con cada atributo en la matriz $X$ y $H$ es la estimación actual.
Los parámetros que se deben entregar a la función de regresión lineal, es la matriz de atributos $X$ y el vector de los valores verdaderos $Y$.
e) Construya una tabla con los pesos y Z-score correspondientes a cada predictor (variable). ¿Qué variables están más correlacionadas con la respuesta? Si usáramos un nivel de significación del 5%. ¿Qué es lo que observa y cuál puede ser la causa?
Los coeficientes corresponden a:
In [5]:
# Obtención de los coeficientes
betas = linreg.coef_
col = list(X)
bT = pd.DataFrame([betas], columns=col)
bT
Out[5]:
Los z-scores de los coeficientes corresponde a:
In [6]:
# Obtención de los z-score
from scipy import stats
zscores = stats.zscore(betas)
zT = pd.DataFrame([zscores], columns=col)
zT
Out[6]:
Lo anterior se puede resumir en la siguiente tabla:
In [7]:
table = bT.append(zT).transpose()
table.columns = ["Coeficiente", "Z-Score"]
table
Out[7]:
Usando un nivel de significancia de $5\%$, tenemos que son relevantes aquellos atributos cuyo $|z|$ sea mayor a las colas de la distribución $\mathcal{N}(\mu, \sigma)$ donde $\mu$ la media de los betas y $\sigma$ la desviación estándar. Como queremos un intervalo de confianza cerrado de igual magnitud, el problema se traduce a calcular:
$2P(Z \leq z) = 0.05$
Tal y como se explica en el libro guía, se puede testear utilizando la distribución $t_{(N - p - 1)}$ donde $N$ es la cantidad de datos utilizados y $p$ la cantidad de atributos, debido a la marginalidad de la diferencia entre la distribución normal y dicha distribución.
Obtenemos entonces:
In [8]:
stats.t.ppf([1-0.025], len(Xtrain)-len(col)-1)
Out[8]:
Lo cual indica que son relevantes aquellas variables con z-score $|z| > 1.96$. Esto es, las variables sqft_living y sqft_above.
Se puede observar que de todos los atributos existentes, la mayoría tiene z-scores marginales en comparaciñon a los atributos sqft_living, sqft_above y sqft_basement, por lo que su aporte es también marginal en la predicción del precio de la vivienda.
f) Proponga un método para corregir lo observado (Hint: inspírese en los métodos de feature engineering de las siguiente secciones). Verifíquelo mediante los Z-score presentados en la pregunta e).
Basandonos en la sección anterior, podemos eliminar todos los atributos con z-score irrelevantes para la regresión lineal. Si consideramos a $z \rel 0.1$ como el umbral de selección, tenemos como resultado:
In [9]:
import sklearn.linear_model as lm
X2 = df_scaled.iloc[:, 1:]
X2 = X2.iloc[:,[2, 9, 10]] #use .ix instead, in older pandas version
print(list(X2))
N = X2.shape[0]
X2.insert(X2.shape[1], 'intercept', np.ones(N))
y2 = df_scaled['price']
#mascara estatica con el 70% de los datos
mascara = np.zeros(len(X2))
limit = int(len(X2)*0.7)
mascara[:limit] = 1
istrain2 = mascara== 1
Xtrain2 = X2[istrain2]
ytrain2 = y2[istrain2]
Xtest2 = X2[np.logical_not(istrain2)]
ytest2 = y2[np.logical_not(istrain2)]
linreg2 = lm.LinearRegression(fit_intercept = False)
linreg2.fit(Xtrain2, ytrain2)
betas2 = linreg2.coef_
col2 = list(X2)
pd.DataFrame([betas2], columns=col2)
zscores2 = stats.zscore(betas2)
pd.DataFrame([zscores2], columns=col2)
Out[9]:
Considerando la distribución $t_{(N - k - 1)}$ donde $k$ son los atributos filtrados. Tenemos entonces:
In [10]:
stats.t.ppf([1-0.025], len(Xtrain2)-len(col2)-1)
Out[10]:
Esto implica que los atributos seleccionados son similarmente importantes en la predicción del precio.
g) Estime el error de predicción del modelo usando validación cruzada con un número de folds igual a K = 5 y K = 10. Recuerde que para que la estimación sea razonable, en cada configuración (fold) deberá reajustar los pesos del modelo. Mida el error real del modelo sobre el conjunto de pruebas, compare y concluya.
Asumiendo que se debe utilizar el modelo con la selección de atributos seleccionada anteriormente, se tiene:
In [11]:
# Using the test sample conaining only selected predictors
yhat_test = linreg.predict(Xtest) # Predict
mse_test = np.mean(np.power(yhat_test - ytest, 2))
# Form the matrix from the RAW set
Xm = Xtrain.as_matrix()
ym = ytrain.as_matrix()
from sklearn.model_selection import KFold
def doKFold(Xt, yt, K):
# Using K = 5
kf = KFold(n_splits=K)
mse_cv = 0
for train_indeces, test_indeces in kf.split(Xm):
# Select only the most relevant attributes
# In the previous section we chose attributes with index 2, 9 and 10
# as the most relevant.
Xtrain_k = Xtrain.iloc[train_indeces].as_matrix()
Xtest_k = Xtrain.iloc[test_indeces].as_matrix()
ytrain_k = ytrain.iloc[train_indeces].as_matrix()
ytest_k = ytrain.iloc[test_indeces].as_matrix()
# Perform linear regression for this fold
linreg = lm.LinearRegression(fit_intercept = False)
linreg.fit(Xtrain_k, ytrain_k)
# Get the yhat for our test set
yhat_val = linreg.predict(Xtest_k)
# Calculate the MSE for the error
mse_fold = np.mean(np.power(yhat_val - ytest_k, 2))
mse_cv += mse_fold
mse_cv = mse_cv / K
return mse_cv
{"mse5": doKFold(Xtrain, ytrain, 5), "mse10": doKFold(Xtrain, ytrain, 10), "mse_test": mse_test}
Out[11]:
Se puede observar que los errores obtenidos para cross validation con $K = 5$ y $K = 10$ son muy similares y tienen un valor aproximado de $0.0647$. El error obtenido del conjunto de prueba sobre los datos sin realización de selección de atributos corresponde a $0.065$ que es prácticamente el mismo.
Esto sucede debido a la gran cantidad de datos disponibles para el entrenamiento. Eligiendo todo el subconjunto de atributos, tenemos que $p$ es al menos 3 ordenes de magnitud menor que $N$, lo que permite tener una gran precisión utilizando cualquier método de validación.
h) Mida los errores de predicción para cada dato de entrenamiento. Utilizando un “quantile-quantile plot” determine si es razonable la hipótesis de normalidad sobre los residuos del modelo.
In [12]:
import scipy.stats as stats
import matplotlib.pyplot as plt
ypredicted = linreg.predict(Xtrain)
error = ypredicted - ytrain.as_matrix()
stats.probplot(error, dist='norm', plot=plt)
plt.ylabel('Error')
plt.title('Quantile-Quantile Plot')
plt.show()
Con el gráfico de Quantile-Quantile se busca analizar las diferencias entre la distribución de una población aleatoria y, en este caso, la distribución normal. Se puede apreciar que los valores están distribuidos en el gráfico de manera aproximadamente lineal con respecto a la recta, lo que es suficiente para respaldar la hipótesis de que los datos están normalmente distribuidos.
Utilizando el dataframe de la actividad anterior,
a) Construya una función que implemente Forward Step-wise Selection (FSS). Es decir, partiendo con un modelo sin predictores (variables), agregue un predictor a la vez, re-ajustando el modelo de regresión en cada paso. Para seleccionar localmente una variable, proponga/implemente un criterio distinto al utilizado en el código de ejemplo. Construya un gráfico que muestre el error de entrenamiento y el error de pruebas como función del número de variables en el modelo. Ordene el eje x de menor a mayor.
Para la implementación de FSS, se ha seleccionado el criterio del promedio del valor absoluto dado por:
$\sum_i^N\frac{|Y - \hat{f}(X)|}{N}$
como una forma alternativa para discriminar a los distintos candidatos.
A continuación se presenta el algoritmo y un gráfico entre el número de variables y el error:
In [13]:
def fss_cuadratic(x, y, names_x, k = 10000):
p = x.shape[1]-1 #p is the total number of attributes
k = min(p, k) #k is the maximum number of parameters to choose
names_x = np.array(names_x) #this is the names of the columns
remaining = list(range(0, p)) #Amount of comparable attributos left
selected = [p] #First, choose the last parameter in the list
current_score = best_new_score = 0.0 #Initialize score
points = []
scores = []
while remaining and len(selected)<=k :
score_candidates = []
for candidate in remaining:
model = lm.LinearRegression(fit_intercept=False) #Init a linear regression model
indexes = selected + [candidate] # Only use the following parameters
x_train = x.iloc[:, indexes]
predictions_train = model.fit(x_train, y).predict(x_train)
residuals_train = predictions_train - y
mse_candidate = np.mean(np.power(residuals_train, 2)) # Promedio del valor absoluto
score_candidates.append((mse_candidate, candidate))
score_candidates.sort()
score_candidates[:] = score_candidates[::-1]
best_new_score, best_candidate = score_candidates.pop()
remaining.remove(best_candidate)
selected.append(best_candidate)
points.append(len(selected))
scores.append(best_new_score)
return selected, points, scores
def fss_abs(x, y, names_x, k = 10000):
p = x.shape[1]-1 #p is the total number of attributes
k = min(p, k) #k is the maximum number of parameters to choose
names_x = np.array(names_x) #this is the names of the columns
remaining = list(range(0, p)) #Amount of comparable attributos left
selected = [p] #First, choose the last parameter in the list
current_score = best_new_score = 0.0 #Initialize score
points = []
scores = []
while remaining and len(selected)<=k :
score_candidates = []
for candidate in remaining:
model = lm.LinearRegression(fit_intercept=False) #Init a linear regression model
indexes = selected + [candidate] # Only use the following parameters
x_train = x.iloc[:, indexes]
predictions_train = model.fit(x_train, y).predict(x_train)
residuals_train = predictions_train - y
mse_candidate = np.mean(np.absolute(residuals_train)) # Promedio del valor absoluto
score_candidates.append((mse_candidate, candidate))
score_candidates.sort()
score_candidates[:] = score_candidates[::-1]
best_new_score, best_candidate = score_candidates.pop()
remaining.remove(best_candidate)
selected.append(best_candidate)
points.append(len(selected))
scores.append(best_new_score)
return selected, points, scores
names_regressors = X.columns[:-1] #without intercept
(selected, points, scores) = fss_abs(X,y,names_regressors)
plt.plot(points, scores)
plt.axis()
plt.title('Error en función al número de predictores')
plt.show()
Como se puede apreciar en el gráfico, a medida que aumenta el número de atributos seleccionados, disminuye el error de la predicción hasta que aproximadamente en $p = 12$ la ganancia de precisión ya no es tan significativa. Esto es relevante al momento de seleccionar atributos cuidando de no caer en overfitting, ya que implica que existe un número ideal o óptimo de atributos a seleccionar.
A continuación se presenta una comparación entre la función cuadratica y absoluta para el error:
In [14]:
(selected, points, scores) = fss_abs(X,y,names_regressors)
plt.plot(points, scores, 'ro')
plt.axis()
plt.title('Absolute error')
plt.show()
(selected, points, scores) = fss_cuadratic(X,y,names_regressors)
plt.plot(points, scores, 'ro')
plt.axis()
plt.title('Cuadratic error')
plt.show()
Se puede notar claramente que utilizando un criterio de error cuadratico se obtienen errores de menor magnitud debido a la naturaleza de la distancia euclidiana. Sin embargo, se puede también notar que ambos siguen la misma tendencia.
Utilizando el dataframe de la actividad anterior,
a) Ajuste un modelo lineal utilizando “Ridge Regression”, es decir, regularizando con la norma $\ell_2$. Utilice valores del parámetro de regularización $\lambda$ en el rango $[10^7, 10^1]$, variando si estima conveniente. Construya un gráfico que muestre los coeficientes obtenidos como función del parámetro de regularización. Describa lo que observa. (WARNING: Note que la línea 3 y el primer argumento en la línea 9 son críticos).
In [15]:
from sklearn.linear_model import Ridge
import matplotlib.pylab as plt
X2 = X.drop('intercept', axis=1,inplace=False)
Xtrain = X2[istrain]
ytrain = y[istrain]
names_regressors = X2.columns
alphas_ = np.logspace(7,1,base=10)
coefs = []
model = Ridge(fit_intercept=True,solver='svd')
for a in alphas_:
model.set_params(alpha=a)
model.fit(Xtrain, ytrain)
coefs.append(model.coef_)
ax = plt.gca()
for y_arr, label in zip(np.squeeze(coefs).T, names_regressors):
plt.plot(alphas_, y_arr, label=label)
plt.legend()
ax.set_xscale('log')
ax.set_xlim(ax.get_xlim()[::-1]) # reverse axis
plt.title('Regularization Path RIDGE')
plt.axis('tight')
plt.legend(loc=2)
plt.show()
Para poder entender el gráfico generado, es necesario entender como funciona Ridge Regression. Este algoritmo penaliza los coeficientes de la regresión lineal a través de un coeficiente denominado $\lambda$. Esto permite "suavizar" la diferencia de valor entre los coeficientes al imponer restricciones sobre su valor máximo.
Cuando se elije un valor de $\lambda$ pequeño, se puede estar en presencia de underfitting, mientras que si es muy pequeño se sufre de overfitting. Por tanto, es importante elegir valores de $\lambda$ adecuados.
Se puede observar en el gráfico que valores de $\lambda$ cercanos a 0 (o bien, cercanos a $10^1$) generan coeficientes muy dispersos con alta varianza y bias. A medida que se aumenta el valor de $\lambda$ los valores de los coeficientes van disminuyendo proporcionalmente de forma que se genera una disminución en la varianza y el bias. Alrededor de un coeficiente de $10^6$, se tiene que los coeficientes comienzan a ser muy similares entre sí mismos, produciendo el efecto de overtfitting que disminuye la capacidad de generalización del modelo entrenado.
Se puede confirmar además que los atributos relevantes corresponden a sqft_lot (Metros cuadrados de la vivienda), grade y bathrooms, debido a que son los que más aportan a la varianza.
b) Ajuste un modelo lineal utilizando el método “Lasso”, es decir, regularizando con la norma $\ell_1$. Utilice valores del parámetro de regularización $\lambda$ en el rango $[10^0, 10^{-3}]$. Para obtener el código, modifique las líneas 7 y 9 del ejemplo anterior. Construya un gráfico que muestre los coeficientes obtenidos como función del parámetro de regularización. Describa lo que observa. ¿Es más efectivo Lasso para seleccionar atributos?
In [16]:
from sklearn.linear_model import Lasso
alphas_ = np.logspace(0,-3,base=10)
X2 = X.drop('intercept', axis=1,inplace=False)
Xtrain = X2[istrain]
ytrain = y[istrain]
names_regressors = X2.columns
coefs = []
model = Lasso(fit_intercept=True)
for a in alphas_:
model.set_params(alpha=a)
model.fit(Xtrain, ytrain)
coefs.append(model.coef_)
ax = plt.gca()
for y_arr, label in zip(np.squeeze(coefs).T, names_regressors):
plt.plot(alphas_, y_arr, label=label)
plt.legend()
ax.set_xscale('log')
ax.set_xlim(ax.get_xlim()[::-1]) # reverse axis
plt.title('Regularization Path LASSO')
plt.axis('tight')
plt.legend(loc=2)
plt.show()
A diferencia de Ride Regression, Lasso permite anular coeficientes para cierto valores de $\lambda$. A esto se le denomina "soft thresholding". Se puede observar que para coeficientes de $\lambda$ pequeños se tiene una gran varianza entre los coeficientes de la regresión lineal. A medida que se aumenta $\lambda$, comienza a disminuir la la varianza entre los coeficientes, pero también la gran mayorìa de ellos comienza a desaparecer. Esto permite reducir la cantidad de atributos requeridos para realizar la estimación, disminuyendo las probabilidades de overfitting si se selecciona un valor de $\lambda$ adecuado.
c) Escogiendo uno de los dos métodos regularizadores anteriores, especificando el porqué, construya un gráfico que muestre el error de entrenamiento y el error de pruebas como función del parámetro de regularización. Discuta lo que observa.
Se utiliza Lasso para aprovechar la oportunidad de reducir ciertos coeficientes a 0 (los que están más cerca del valor $\lambda = 10^{-3}$) en lugar del suavizado que realiza Ridge Regression.
In [17]:
Xtest = X2[np.logical_not(istrain)]
ytest = y[np.logical_not(istrain)]
alphas_ = np.logspace(0,-3,base=10)
coefs = []
model = Lasso(fit_intercept=True)
mse_test = []
mse_train = []
for a in alphas_:
model.set_params(alpha=a)
model.fit(Xtrain, ytrain)
yhat_train = model.predict(Xtrain)
yhat_test = model.predict(Xtest)
mse_train.append(np.mean(np.power(yhat_train - ytrain, 2)))
mse_test.append(np.mean(np.power(yhat_test - ytest, 2)))
ax = plt.gca()
ax.plot(alphas_,mse_train,label='train error ridge')
ax.plot(alphas_,mse_test,label='test error ridge')
plt.legend(loc=1)
ax.set_xscale('log')
ax.set_xlim(ax.get_xlim()[::-1])
plt.show()
Se puede notar que a medida que se disminuye el valor del parámetro $\lambda$, comienzan a aumentar tanto el error de entrenamiento como el error de prueba. Esto ocurre debido a que en Lasso muchos de los coeficientes se anulan, evitando que sus atributos asociados aporten información a la predicción. Entre $10^0$ y $10^{-1}$ se llega a un punto crítico en donde se obtiene un error máximo debido a la elminicación de la mayoría de los atributos.
d) Estime el valor del parámetro de regularización en alguno de los modelos anteriores haciendo uso de la técnica validación cruzada.
Utilizando Lasso, se tiene que:
In [18]:
def MSE(y,yhat): return np.mean(np.power(y-yhat,2))
Xm = Xtrain.as_matrix()
ym = ytrain.as_matrix()
from sklearn.model_selection import KFold
kf = KFold(n_splits=10)
best_cv_mse = float("inf")
alphas_ = np.logspace(0,-3,base=10)
model = Lasso(fit_intercept=True)
for a in alphas_:
model.set_params(alpha=a)
mse_list_k10 = [MSE(model.fit(Xm[train], ym[train]).predict(Xm[val]), ym[val]) for train, val in kf.split(Xm)]
if np.mean(mse_list_k10) < best_cv_mse:
best_cv_mse = np.mean(mse_list_k10)
best_alpha = a
print("BEST PARAMETER=%f, MSE(CV)=%f"%(best_alpha,best_cv_mse))
Por lo tanto, el parámetro ideal para Lasso corresponde a $10^{-3}$ que corresponde al intercepto entre el error de entrenamiento y el error de testing de menor valor, coincidiendo con el gráfico anterior.
En esta sección se presentarán dos muestras del dataframe utilizado en la actividades anteriores, donde cada una de estas tiene una propiedad distinta ya que son muestreadas en función del valor a predecir (logaritmo del precio de la casa). Por una parte se tiene una pequeña muestra A, la cual es extraída directamente de los datos con los que se trabaja (manteniendo la distribución de esta) y la muestra B, es generada con el propósito de que en cada intervalo del rango de valores haya la misma cantidad de datos aproximadamente (simulando una distribución uniforme). El objetivo es familiarizarse con el concepto de Transfer Learning.
En el siguiente código se generan las dos muestras con las que se trabajará.
In [19]:
import pandas as pd
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn import linear_model
df = pd.read_csv("kc_house_data.csv")
df.drop(['id','date','zipcode',],axis=1,inplace=True)
scaler = StandardScaler()
df_scaled = pd.DataFrame(scaler.fit_transform(df), columns=df.columns)
df_scaled['price'] = np.log(df['price'])
df_A = df_scaled.sample(1000,random_state=11)
frames = []
valor = df_scaled.price
length = 0.3
for z in np.arange(int(np.min(valor)),int(np.max(valor))+1,length):
#un maximo de 100 datos por intervalo
aux = df_scaled[(df_scaled.price >= z) & (df_scaled.price < z+length)].head(100)
frames.append(aux)
df_B = pd.concat(frames).sample(1000,random_state=11) #crea el dataframe
a) Cree el conjunto de entrenamiento y otro de validación para trabajar cada muestra mediante la técnica de hold out validation.
In [20]:
X_A = df_A.iloc[:,1:].values
y_A = df_A.price
X_B = df_B.iloc[:,1:].values
y_B = df_B.price
from sklearn.model_selection import train_test_split
Xtrain_A,Xval_A,ytrain_A,yval_A = train_test_split(X_A, y_A, test_size=0.3, random_state=42)
Xtrain_B,Xval_B,ytrain_B,yval_B = train_test_split(X_B, y_B, test_size=0.3, random_state=42)
b) Evalúe los dos modelo de regresión lineal que se generan al entrenar con cada muestra. Mida el error de cada modelo sobre ambos conjuntos de validación (A y B). Explique lo que observa.
In [21]:
lm_A = linear_model.LinearRegression(fit_intercept=False)
model_A = lm_A.fit(Xtrain_A, ytrain_A)
predictions_A_test_A = model_A.predict(Xval_A)
predictions_A_test_B = model_A.predict(Xval_B)
lm_B = linear_model.LinearRegression(fit_intercept=False)
model_B = lm_B.fit(Xtrain_B, ytrain_B)
predictions_B_test_A = model_B.predict(Xval_A)
predictions_B_test_B = model_B.predict(Xval_B)
# print(predictions_A, predictions_B)
msa_a_test_a = np.mean(np.power(predictions_A_test_A - yval_A, 2))
msa_a_test_b = np.mean(np.power(predictions_A_test_B - yval_B, 2))
msa_b_test_a = np.mean(np.power(predictions_B_test_A - yval_A, 2))
msa_b_test_b = np.mean(np.power(predictions_B_test_B - yval_B, 2))
score_a_test_a = model_A.score(Xval_A, yval_A)
score_a_test_b = model_A.score(Xval_B, yval_B)
score_b_test_a = model_B.score(Xval_A, yval_A)
score_b_test_b = model_B.score(Xval_B, yval_B)
print("MSA Train A Test A:", msa_a_test_a)
print("MSA Train A Test B:", msa_a_test_b)
print("MSA Train B Test A:", msa_b_test_a)
print("MSA Train B Test B:", msa_b_test_b)
print("MSA Average Train A", (msa_a_test_a+msa_a_test_b)/2)
print("MSA Average Train B", (msa_b_test_a+msa_b_test_b)/2)
print("-----------------")
print("Score Train A Test A:", score_a_test_a)
print("Score Train A Test B:", score_a_test_b)
print("Score Train B Test A:", score_b_test_a)
print("Score Train B Test B:", score_b_test_b)
print("Score Average Train A", (score_a_test_a+score_a_test_b)/2)
print("Score Average Train B", (score_b_test_a+score_b_test_b)/2).
Se puede observar por simple inspección que se tiene en promedio un error cuadratico menor para la muestra B. Además, se observa que al utilizar un conjunto de prueba distinto al asociado con cada modelo (es decir, usar el conjunto de testing de B en el modelo A y viceversa) produce errores mayores en comparación a utilizar el conjunto de prueba original. Esto se debe a que inveitablemente cada modelo tenderá a producir overfitting para la distribución especifica en que estén los datos de entrnamiento. Si los datos de testing siguen la misma distribución, entonces tendrán menor error.
c) Si tuviera que elegir uno de los dos modelos anteriores para trabajar con data futura, ¿Cuál eligiría y por qué?
No es tan facil realizar una decisión sin saber a priori la distribución real de los datos. Por ejemplo, tenemos que el modelo B tiene un error mucho menor en conjuntos de testing con su misma distribución, pero arroja un error mucho mayor (posee más overfitting) para conjuntos con una distribución distinta. El modelo A por otro lado, tiene mayor error para el conjunto de testing de su misma distribución, pero arroja un error menor para conjuntos con distribuciones distintas.
Aún así, y considerando que en la práctica obtener conjuntos uniformemente distribuidos de precios de casas no es razonable, consideramos que es mejor utilizar el modelo A para reducir la perdida de generalización.
En el área de la salud, diagnosticar a una persona de una enfermedad de forma rápida y correcta puede llegar a salvarle la vida. Los encargados de realizar estos diagnósticos, son médicos que, observando exámenes y ciertos indicadores, pueden concluir qué enfermedad presenta el paciente. Si el medico se llegase a equivocar, aparte de que el paciente pueda perder la vida, el medico podría ser demandado por negligencia arriesgando años de cárcel o pagar sumas de dinero considerable, es por estas razones que es importante no cometer errores. Pongámonos en el contexto de que usted es contratado para generar un modelo que prediga si es que un paciente presenta una enfermedad cardiaca a partir de ciertos indicadores, tales como la edad, sexo, presión sanguínea, nivel de glicemia, etc.
Como ayuda se le indica que la variable de máximo ritmo cardíaco alcanzado (maximum heart rate achieved) es un buen indicador de detección de enfermedades cardíacas. Por lo que el objetivo es predecir el comportamiento de esta variable en función de las otras, y con esto detectar qué tan distante es el valor real al valor predecido para así luego detectar los posibles outliers (enfermos), que en sí corresponden a pacientes que tienen un comportamiento anormal al resto.
a) Lea el archivo de datos, cárguelos en un dataframe o matriz, luego divida el dataframe en dos, un dataframe de entrenamiento (70% Datos) y un dataframe de prueba (30% Datos).
In [22]:
import pandas as pd
from sklearn import linear_model
from matplotlib import pyplot as plt
from sklearn.model_selection import train_test_split
headers = ['age','sex','chest_pain','blood_p','serum','blood_s','electro','max_heart','angina','oldpeak','slope','vessel','thal','normal']
df = pd.read_csv('heart.csv', header=None, names=headers, sep=',')
y = df['max_heart']
# train_n = int(len(df)*0.7)
# X_train, X_test = df[-train_n:], df[:-train_n]
# y_train, y_test = y[-train_n:], y[:-train_n]
X_train, X_test, y_train, y_test = train_test_split(df[headers], y, test_size=0.3, random_state=42)
X_train_normal = X_train['normal']
X_train_heart = X_train['max_heart']
X_train_class = X_train[['normal', 'max_heart']]
X_train.drop(['max_heart','normal'],axis=1,inplace=True)
X_test_normal = X_test['normal']
X_test_heart = X_test['max_heart']
X_test.drop(['max_heart','normal'],axis=1,inplace=True)
X_train
Out[22]:
b) Realice una regresión lineal y defina usted una frontera de decisión (umbral) para poder determinar si es que estamos en presencia o no de una enfermedad cardíaca. Mida su desempeño con ambos conjuntos de datos.
Para obtener un umbral relativamente aceptable, utilizaremos el promedio de valores anormales menos la mitad de la desviación estándar de los valores. Se tiene entonces:
In [23]:
from sklearn.metrics import accuracy_score
lm = linear_model.LinearRegression()
model = lm.fit(X_train, y_train)
predictions_train = model.predict(X_train)
predictions_test = model.predict(X_test)
# Calculo del umbral
sick = X_train_class[X_train_class['normal'] == 1]['max_heart']
limit = sick.mean()-sick.std()/2
print("Umbral: ", limit)
compare_train = pd.DataFrame({'Y': y_train, 'y_': predictions_train})
compare_test = pd.DataFrame({'Y': y_test, 'y_': predictions_test})
mse_train = np.mean(np.power(predictions_train - y_train, 2))
mse_test = np.mean(np.power(predictions_test - y_test, 2))
print("mse_train:", mse_train)
print("mse_test", mse_test)
# Clasificación original
y_train_outlier = X_train_normal
y_test_outlier = X_test_normal
def heart_to_class(val):
if(val >= limit):
return 1
return 2
# Aplicar clasificación según umbral
y_train_predict_outlier = compare_train['y_'].apply(heart_to_class)
y_test_predict_outlier = compare_test['y_'].apply(heart_to_class)
print("Score: ", accuracy_score(y_train_outlier,y_train_predict_outlier))
print("Score: ", accuracy_score(y_test_outlier,y_test_predict_outlier))
En donde se tiene un puntaje de 80% para el conjunto de entrenamiento y 71% para el de testing.
A continuación se prueba con distintos umbrales:
In [24]:
def attemptClassification(limit):
def heart_to_class(val):
if(val >= limit):
return 1
return 2
# Aplicar clasificación según umbral
y_train_predict_outlier = compare_train['y_'].apply(heart_to_class)
y_test_predict_outlier = compare_test['y_'].apply(heart_to_class)
return(accuracy_score(y_train_outlier,y_train_predict_outlier), accuracy_score(y_test_outlier,y_test_predict_outlier))
print("2std", attemptClassification(sick.mean()-sick.std()*2))
print("std", attemptClassification(sick.mean()-sick.std()))
print("std/2", attemptClassification(sick.mean()-sick.std()/2))
print("std/4", attemptClassification(sick.mean()-sick.std()/4))
print("std/8", attemptClassification(sick.mean()-sick.std()/8))
print("std/16", attemptClassification(sick.mean()-sick.std()/16))
std_coef = pd.DataFrame({'coef': [0.1, 0.2, 0.3, 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 16, 32] })
std_coef['coef'] = std_coef['coef'].apply(lambda x: attemptClassification(sick.mean()-sick.std()*x)[0])
std_coef_test = pd.DataFrame({'coef': [0.1, 0.2, 0.3, 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 16, 32] })
std_coef_test['coef'] = std_coef_test['coef'].apply(lambda x: attemptClassification(sick.mean()-sick.std()*x)[1])
plt.plot([0.1, 0.2, 0.3, 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 16, 32], std_coef['coef'].as_matrix(), 'ro')
plt.axis()
plt.title('Train score')
plt.show()
plt.plot([0.1, 0.2, 0.3, 0.5, 1, 2, 3, 4, 5, 6, 7, 8, 16, 32], std_coef_test['coef'].as_matrix(), 'ro')
plt.axis()
plt.title('Test score')
plt.show()
In [ ]:
In [ ]: