<jupyter_start><jupyter_text>Cargamos las librerías iniciales<jupyter_code>import pandas as pd
import matplotlib.pyplot as plt 
import seaborn as sns
import numpy as np

%matplotlib inline<jupyter_output><empty_output><jupyter_text>Descripción de ColumnasLas columnas del dataset describen aspectos físicos del asteroide como también de su orbita con respecto a la Tierra. Todas las características de órbita son con respecto a la Tierra.* Semi-Major Axis (a): Se refiere a la mitad del diametro más largo de la elipse* Eccentricity (e): Se refiere a la cantidad en la que la orbita del asteroide se desvía de la forma de un círculo* Inclination with respect to x-y ecliptic plane (deg): "Although all the planets and asteroids follow elliptical orbits around the Sun (obeying Kepler’s First Law), these orbits do not all lie in the same plane – they are usually tilted with respect to each other. As Earth-bound humans, we have adopted the plane in which the Earth moves around the Sun (the ecliptic) as our reference plane for the Solar System. With this convention, the Earth has an orbital inclination of zero degrees, and the orbital inclinations of other Solar System bodies are measured relative to this (for example, Mars has an orbital inclination of 1.85o, Mercury: 7.00o and Pluto: 17.15o)." Tomado de https://astronomy.swin.edu.au/cosmos/O/orbital+inclination. En nuestro caso sería una referencia relativa al asteroide.* Longitude of the Ascending Node (om): "The point where the orbit passes through the reference plane in an “upward” direction is called the ascending node. Thus, the second orientation angle of the orbit is the angle between the ascending node and some reference direction in the reference plane, as seen from the Sun. In our solar system, the reference direction is usually taken to be the direction toward the vernal equinox. This angle is known as the longitude of the ascending node, Ω." Tomado de: https://www.sciencedirect.com/topics/physics-and-astronomy/longitude-of-the-ascending-node* Argument of Perihelion (w): Es la distancia angular del ascending node hasta la el perhilion medido en el plano orbital.* Periphelion Distance (q): Es la distancia de la orbita del asteroide en la cual se encuentra lo más cercano al sol.* Aphelion Distance (ad): Es la distancia de la orbita del asteroide en la cual se encuentra lo más alejada del sol.* Orbital Period in Years (per_y): Se refiere al tiempo en que se tarda el asteroide en completar su orbita alrededor de la Tierra en años.* Data Arc-Span (data_arc): Sin descripción* Orbit Condition Code (condition_code): Se refiere a un parámetro de incertidumbre sobre la órbita del asteroide* Number of observation used (n_obs_used): Número de observaciones que se le hicieron al asteroide para establecer las demás características.* Absolute Magnitude Parameter (H): Se refiere a la magnitud en la cual una estrella aparecería estar si se situara a 10 parsecs de un asteroide (1 parsec = 30 trillones de km aprox = 3 años luz)* Near Earth Object (neo): Es un objeto celestial que pasa relativamente cerca a la órbita de la Tierra* Physically Hazardous Asteroid (pha) VARIABLE OBJETIVO: Es un asteroide que hace trayectos cercanos a la órbita de la Tierra y tiene el tamaño suficiente como para poder causar daños severos en caso de impacto con nuestro planeta.* Diameter of Asteroid in km (diameter): Diámetro del asteroide en kilómetros* Object bi/tri axial ellipsoid dimensions km (extent): Sin descripción* Geometric Albedo (albedo): "...the ratio of its total brightness at zero phase angle to that of an idealised fully reflecting, diffusively scattering (Lambertian) disk with the same cross-section." Tomado de https://space.fandom.com/wiki/Geometric_albedo* Rotation Period in hours (rot_per): Tiempo en horas que le tarda al asteroide hacer una revolución en su propio eje.* Standard Gravitational Parameter (GM): Producto de la masa y la constante gravitacional* Color Index B-V magnitude difference (BV): Se refiere a una escala de medición de temperatura del cuerpo celestial. Números más pequeños indican que es caliente, mientras que números más grandes indican que es más caliente.* Color Index U-B magnitude difference (UB): Se refiere a una escala de medición de temperatura del cuerpo celestial. Números más pequeños indican que es caliente, mientras que números más grandes indican que es más caliente.* Color Index I-R magnitude difference (IR): Se refiere a una escala de medición de temperatura del cuerpo celestial. Números más pequeños indican que es caliente, mientras que números más grandes indican que es más caliente.* Spectral Taxonomic type SMASSII (spec_B): Una clasificación dada al color que emite el asteroide* Spectral Taxonomic type Tholen (spec_T): Una clasificación dada al color que emite el asteroide* Magnitude Slope Parameter (G): No hay descripción* Earth Minimum Orbit Intersection Distance au (moid): "The minimum orbital intersection distance (MOID) is used as a measure to assess potential close approaches and collision risks between astronomical objects" Tomado de https://www.aanda.org/articles/aa/full_html/2020/01/aa36502-19/aa36502-19.html* Asteroid Orbit Class (class): MBA se refiere a Main Belt Asteroid, son asteroides que se encuentran en el anillo de asteroides que se encuentra entre Marte y Júpiter. OMB se refiere a Outer Main Belt Asteroid, son asteroides que se encuentran en el anillo de asteroides que se encuentran hasta Neptuno. Existen varias diferentes clasificaciones.* Mean motion deg/d (n): No hay descripción* Orbital Period d (per): Días que se tarda el asteroide en completar su orbita alrededor de la Tierra* Mean Anomaly deg (ma): "The mean anomaly is the angle between lines drawn from the Sun to the perihelion B and to a point (not shown) moving in the orbit at a uniform rate corresponding to the period of revolution of the planet." Tomado de https://www.britannica.com/science/anomaly-astronomyref105659<jupyter_code>df = pd.read_csv('Asteroid_Updated.csv')
df.head()
df.shape
df.dtypes<jupyter_output><empty_output><jupyter_text>Diameter debe ser una columna númerica por lo que haremos el cambio de data type correspondiente. No está incluído el error que salió previo a la elaboración final del notebook pero ocurrió un error al tratar de convertir la columna a números ya que era de texto y en algunos campos tenían espacios. Por lo mismo, previo a convertirla númerica le quité todos los espacios en blanco de las observaciones para que la función logrará hacer la transformación.<jupyter_code>df['diameter'] = df['diameter'].str.strip()
df["diameter"] = pd.to_numeric(df["diameter"], downcast="float")
df.shape
df['pha'].value_counts()<jupyter_output><empty_output><jupyter_text>Explicación de distribución de variable objetivo PHA (Physically Hazardous Asteroid)Para tener un poco de contexto, actualmente el número de asteroides, meteoritos y cometas que se consideran peligrosos para la humanidad son un número bastante bajo. Existen innumerables objetos del mismo tipo en nuestro sistema solar. Residuos de colisiones de planetas dentro de nuestro sistema solar, como también de objetos fuera de nuestra vecindad local. Teniendo esto en cuenta, la NASA ha declarado públicamente que el número de objetos que son peligrosos para la Tierra son pocos. El asteroide más notorio actualmente es Apophis, tiene una longitud/diametro de 340 metros y se había sugerido que iba a tener una trayectoria sumamente cercana a la Tierra en el año 2029. Tras varios análisis descartaron la posibilidad y calcularon que en el año 2068 tiene una probabilidad mínima de poder impactar con la Tierra. Con esto en mente, se puede justificar el bajo número de asteroides considerados como potencialmente peligrosos. El tamaño seguramente dicta lo suficiente para considerarse peligroso, pero se debe de tomar en cuenta características de su trayectoria ya que incluso asteroides de a penas cientos de metros son capaces de poder destruir ciudades enteras si se encuentran en la trayectoria justa. El tamaño necesario para erradicar la vida en la Tierra se ha estimado que tiene que ser de aproximadamente 96 km.  Quitamos la columna pha puesto a que vamos a utilizar nuestra propia variable objetivo, ya que sería redundante tener dos columnas que describan algo similar<jupyter_code>df.drop('pha', inplace=True, axis=1)<jupyter_output><empty_output><jupyter_text>Redifinir la variable objetivoSi utilizamos a PHA como la variable objetivo, nuestro modelo seguramente va a hacer underperform drásticamente. Por lo mismo, vamos a cambiar nuestra clasificación a identificar todos los asteroides que tienen la capacidad de extinguir a la raza humana. El diametro necesario como para tener implicaciones globales es de 8km. No necesariamente destruiría físicamente toda la Tierra. Lo que ocurriría con este tipo de asteroide es que en el lugar de impacto obviamente dejaría todo devastado, pero por el impacto causaría cambios climáticos como también la difusión de gases tóxicos que se generaron por el impacto. Este tipo de químicos estando presentes en la atmósfera tienen la capacidad de extenderse a todo el mundo, y terminaría matando a todo ser vivo que estuviera en la superficie. Adicionalmente, también mataría la vida marina puesto a que al momento de llover, las nubes estarían desatando todos los químicos que se encontraran en el aire. Estaremos considerando todos los asteroides que superan los 8 km de diámetro como extinguidores, mientras que los que midan menos los consideraremos no catastróficos.  Importancia de crear modelosLos modelos de machine learning tienen que ser capaces de poder predecir qué asteroides son potencialmente peligrosos para la Tierra. Estos modelos van a tener que trabajar con este tipo de data, donde la variable objetivo no va a tener una distribución adecuada y probablemente encuentren más asteroides que no son una amenaza para nosotros. Actualmente no contamos con los sistemas necesarios para actuar y responder a estas amenazas. Sin embargo, la NASA ya está en proceso de crear un sistema de respuesta en el caso un asteroide tenga alta probabilidad de impacto con la Tierra. Su método de respuesta involucra envíar una cápsula espacial considerablemente grande, para impactar el asteroide y tratar de desviar su trayectoria. El asteroide más peligroso no es el que se detecta, son todos los que no se logran detectar y tienen alta probabilidad de impactarnos.  Análisis de NA'sPodemos observar que existen varias columnas con NA's. No sería adecuado aplicar imputación ya que los NA's representan más del 50% de las observaciones. Si se aplicara imputación estaríamos cometiendo un error grave dado a que aunque utilizemos modelos de imputación, el algoritmo estaría prácticamente adivinando o teniendo predicciones poco precisas. La columna con mayor importancia es diámetro ya que es nuestra variable objetivo. No podemos tener NA's en la variable que estamos tratando de predecir. Vamos a poder usar solo el 3% total del dataset.<jupyter_code>df.isnull().mean() * 100<jupyter_output><empty_output><jupyter_text>Necesitamos ponerle énfasis a la variable objetivo, por lo cual hacemos un filtro para obtener todos los registros que no contengan NA's en esa columna.<jupyter_code>df2 = df[df['diameter'].notna()]<jupyter_output><empty_output><jupyter_text>Queremos eliminar las columnas que contengan más del 30% del total de datos de NA's<jupyter_code>df3 = df2.loc[:, df2.isin([' ','NULL',0, np.nan]).mean() < .30]
df3.isnull().sum()
df3.shape<jupyter_output><empty_output><jupyter_text>Categorizar los asteroidesAplicamos primero un filtro para convertir en 1 a todos los asteroides que sean de 8km o más largos de diametro. Después simplemente ponemos que el resto de los datos se vuelvan 0's, identificandolos como no catastróficos.<jupyter_code>df3.loc[df3['diameter'] >= 8, 'diameter'] = 1
df3.loc[df3['diameter'] != 1, 'diameter'] = 0
df3['diameter'].value_counts()<jupyter_output><empty_output><jupyter_text>Podemos ver que la variable objetivo ahora tiene una mejor distribución, a diferencia de la que estabamos utilizando originalmente de *pha*.  Correlation MatrixPodemos visualizar desde ahorita que no existe una correlación tan fuerte de los datos con respecto a la variable objetivo. Las variables que tienen más correlación serían **moid y H.**<jupyter_code>corr = df3.corr()
sns.heatmap(corr, xticklabels = df3.columns, yticklabels=df3.columns, cmap="mako")<jupyter_output><empty_output><jupyter_text>Distribución de los Datos<jupyter_code>df3.hist(figsize=(15,17))<jupyter_output><empty_output><jupyter_text>Normalization FormulaLa mayoría de columnas tienen una distribución no normal, por lo que nos indica que será más adecuado normalizarlos que utilizar un Standard Scaler. (El standard scaler funciona mejor con datos con distribuciones normales). Abajo podrán visualizar cómo funciona cada uno. Normalizer también se le conoce como un MinMaxScaler.      Standardization FormulaEste método hace un rescale a la data para que tenga una media de 0, y una desviación estándar de 1.      Imputation MethodLas imputaciones que haremos para los NA's que todavían están presentes serán hechas por la media, ya que si utilizamos el promedio estaríamos afectando los datos por los outliers que están presentes en la columna de albedo.      Filtrado por QuantilesPara poder asegurar lo más posible un buen set de datos, estaré filtrando por cuartiles la data. Quiero poder eliminar outliers sin tener que perjudicar la cantidad de datos que vamos a trabajar con.<jupyter_code>cols = ['a', 'e', 'i', 'q', 'ad', 'per_y', 'data_arc', 'n_obs_used', 'albedo', 'moid', 'n', 'per']
Q1 = df.quantile(0.05)
Q3 = df.quantile(0.95)
IQR = Q3 - Q1
df4 = df3[~((df3[cols] < (Q1 - 1.5 * IQR)) |(df3[cols] > (Q3 + 1.5 * IQR))).any(axis=1)]
df4.shape
df4['diameter'].value_counts()<jupyter_output><empty_output><jupyter_text>El código de arriba elimina a todos los datos que se encuentren 1.5 veces afuera del intervalo intercuartil. Normalmente los cuartiles serían 25% y 75% pero lo más probable es que nos quedaríamos con un dataset sumamente pequeño. Necesitamos la mayor cantidad de datos posible para poder tener diversidad en las observaciones como también lo suficiente para que el modelo pueda aprender adecuadamente.  Nueva distribución de datos<jupyter_code>df4.hist(figsize=(15,17))<jupyter_output><empty_output><jupyter_text>Logramos mejorar un poco la distribución y ahora se parecen más a una distribución normal. En vez de utilizar normalizer, vamos ahora a utilizar el Standard Scaler para transformar los datos.<jupyter_code>df4.isnull().sum()<jupyter_output><empty_output><jupyter_text>Exportar archivo para reproducibilidadExportamos el archivo reducido para poder utilizarse más adelante en los análisis posteriores. Este lo podrán encontrar adjunto con el entregable sin necesidad de correr los códigos anteriores.<jupyter_code>df4.to_csv('asteroid_reduced.csv', index=False)
df5 = pd.read_csv('asteroid_reduced.csv')<jupyter_output><empty_output><jupyter_text>Separación variables X y Y<jupyter_code>y = df5['diameter']
X = df5.drop('diameter', axis=1)<jupyter_output><empty_output><jupyter_text>Pre procesamiento columnas X's  Pipelines para preprocessingNecesitamos transformar la data para poder entregar al modelo un set de datos que pueda procesar. Utilizamos pipeline poder aplicar one hot encoding a las categóricas (hace dummify las variables para poder representarse en números y que el modelo de machine learning logre procesarlo), median imputer para la columna númerica que tiene NA's, y finalmente un Standard Scaler para rescale los valores<jupyter_code>from sklearn.pipeline import Pipeline
from sklearn.preprocessing import OneHotEncoder
from sklearn.preprocessing import Normalizer
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler


cat_pipeline = Pipeline([
        ('oh_enc', OneHotEncoder())
    ])

num_pipeline = Pipeline([
        ('imputer', SimpleImputer(strategy="median")),
        ('normalize_data', StandardScaler())
    ])

num_col = X.select_dtypes(include='number').columns
cat_col = X.select_dtypes(exclude='number').columns
from sklearn.compose import ColumnTransformer
full_pipeline = ColumnTransformer([
        ("cat", cat_pipeline, cat_col),
        ("num", num_pipeline, num_col)
    ])<jupyter_output><empty_output><jupyter_text>Se juntan los dos transformers para tener el full pipeline necesario para procesar el X_train y X_test  Splitting DataSeparamos el dataset en 70% training y 30% testing, ponemos un random state para reproducibilidad y ponemos el parámetro shuffle como *True*  para que haga shuffle la data<jupyter_code>from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=20, shuffle=True)
X_train_processed = full_pipeline.fit_transform(X_train)
print(X_train_processed.shape)
X_test_processed = full_pipeline.fit_transform(X_test)
print(X_test_processed.shape)<jupyter_output>(6261, 24)<jupyter_text>F1 ScoringEl F1 score es una combinación de las métricas de precision y recall. Nos da una accuracy del modelo.* Precision: Son los positivos verdaderos divididos por la suma de todos los positivos verdaders y todos los negativos verdaderos* Recall: Se refiere a la predicción correcta de los positivos, es la proporción entre los positivos verdaderos y todos los positivos.Si queremos que el modelo tenga más precision, estamos diciendo que preferimos tener errores de un negativo falso. En nuestro caso, sería decir que preferimos cometer el error de clasificar a un asteroide como no peligroso cuando realmente lo es.Si queremos que el modelo tenga más recall, estamos diciendo que preferimos tener errores de un falso positivo. En nuestro caso, sería decir que preferimos clasificar a un asteroide como potencialmente peligroso, cuando realmente no lo es.Ambas métricas son importantes para nosotros, pues es sumamente peligroso no tener en cuenta que un asteroide es peligroso. Sin embargo, tampoco podemos gastar recursos y mandar una nave/cohete hacia un asteroide que no es peligroso ya que estaríamos gastando recursos, esfuerzos y dejándonos vulnerable al impacto de otro asteroide que probablemente sí llegue a la Tierra<jupyter_code>from keras import backend as K
def recall(y_true, y_pred):
    y_true = K.ones_like(y_true) 
    true_positives = K.sum(K.round(K.clip(y_true * y_pred, 0, 1)))
    all_positives = K.sum(K.round(K.clip(y_true, 0, 1)))
    
    recall = true_positives / (all_positives + K.epsilon())
    return recall

def precision(y_true, y_pred):
    y_true = K.ones_like(y_true) 
    true_positives = K.sum(K.round(K.clip(y_true * y_pred, 0, 1)))
    
    predicted_positives = K.sum(K.round(K.clip(y_pred, 0, 1)))
    precision = true_positives / (predicted_positives + K.epsilon())
    return precision

def f1_score(y_true, y_pred):
    precision_2 = precision(y_true, y_pred)
    recall_2 = recall(y_true, y_pred)
    return 2*((precision_2*recall_2)/(precision_2+recall_2+K.epsilon()))<jupyter_output><empty_output><jupyter_text>Modelo Redes NeuronalesNuestra red neuronal tiene 4 capas, de las cuales 3 son hidden layers que procesan y la última es el output layer. El primer layer es de 24 neuronas, ya que ese es la dimensión que contamos con nuestros datos X's. Después de manera consecutiva, se reduce el número de capas por la mitad hasta finalmente tener el output layer. Los hidden layers todos tienen activación relu que es lo estándar al usar redes. La parte más importante es la del output layer. Dependiendo del problema que estemos tratando de solucionar, se debe de escoger la función de activación apropiada para la output layer. En nuestro caso, como estamos prediciendo valores entre 0 y 1, utilizamos la función sigmoide. Abajo podrán ver que traté de utilizar un GridSearchCV para poder encontrar los parámetros óptimos del modelo. Sin embargo, estuve teniendo problemas con el procesamiento como también los resultados. Decidí mejor utilizar los modelos por separado.<jupyter_code>'''
from keras.wrappers.scikit_learn import KerasClassifier
from sklearn.model_selection import GridSearchCV

def create_model(optimizer='adam'):
    model = Sequential()
    model.add(Dense(24, input_dim=24, activation='relu'))
    model.add(Dense(12, activation='relu'))
    model.add(Dense(6, activation='relu'))
    model.add(Dense(1, activation='sigmoid'))
    model.compile(loss='binary_crossentropy', metrics=[f1_score])
    return model 

optimizer = ['SGD', 'adam']
batch_size = [10, 20, 40, 50]
epochs = [25, 50]
param_grid = dict(optimizer=optimizer, batch_size=batch_size, epochs=epochs)
model = KerasClassifier(build_fn=create_model, verbose=2)
grid = GridSearchCV(estimator=model, param_grid=param_grid, n_jobs=-1, cv=3)
grid_result = grid.fit(X_train_processed, y_train)
'''<jupyter_output><empty_output><jupyter_text>Parámetros Loss FunctionDado a que estamos trabajando simplemente con 1's y 0's, utilizamos la loss function de binary crossentropy. Este método compara cada una de las probabilidades predecidas a una clase (1 o 0). Dependiendo de la distancia que tiene con el score real, es como penaliza el modelo.      Optimizador AdamActualiza los pesos de la redes neuronales iterativamente, y mantiene un promedio exponencial que decae de los gradientes anteriores. Este algoritmo se utiliza como un sustito al SGD.<jupyter_code>from keras.models import Sequential
from keras.layers import Dense
from keras import optimizers
import tensorflow as tf
model = Sequential()
model.add(Dense(24, input_dim=24, activation='relu'))
model.add(Dense(12, activation='relu'))
model.add(Dense(6, activation='relu'))
model.add(Dense(1, activation='sigmoid'))
opt = optimizers.Adam(learning_rate=0.01)
model.compile(loss=['binary_crossentropy'], metrics=[f1_score], optimizer=opt)
model.fit(X_train_processed, y_train, epochs=50, batch_size=10, use_multiprocessing=True)
results = model.evaluate(X_test_processed, y_test, batch_size=10, use_multiprocessing=True)
print("test loss, test acc:", results)<jupyter_output>627/627 [==============================] - 0s 506us/step - loss: 0.0954 - f1_score: 0.4584
test loss, test acc: [0.09539930522441864, 0.4583682417869568]<jupyter_text>Podemos ver que nuestro modelo se mantiene constante entre el training y el testing set con un F1 score de 45. Estaremos evaluando el modelo de SDG para ver si podemos mejorar la métrica.  Stochastic Gradient DescentElabora una actualización de parámetros para cada observación del training set y su respectiva *y*. Siempre se inicializan con parámetros random the theta, y después el modelo va iterativamente trabajando en la pendiente de la función hasta que alcanza el punto más bajo posible. Tomado de: https://towardsdatascience.com/stochastic-gradient-descent-clearly-explained-53d239905d31<jupyter_code>model = Sequential()
model.add(Dense(24, input_dim=24, activation='relu'))
model.add(Dense(12, activation='relu'))
model.add(Dense(6, activation='relu'))
model.add(Dense(1, activation='sigmoid'))
opt = tf.keras.optimizers.SGD(learning_rate=0.01)
model.compile(loss='binary_crossentropy', metrics=[f1_score], optimizer=opt)
model.fit(X_train_processed, y_train, epochs=50, batch_size=10, use_multiprocessing=True)
results = model.evaluate(X_test_processed, y_test, batch_size=10, use_multiprocessing=True)
print("test loss, test acc:", results)<jupyter_output>627/627 [==============================] - 0s 525us/step - loss: 0.0928 - f1_score: 0.4600
test loss, test acc: [0.0928306058049202, 0.4600352942943573]
{'fraction_non_alphanumeric': 0.06410581473968897, 'fraction_numerical': 0.01533975659229209, 'mean_word_length': 5.612182173791562, 'pattern_counts': {'":': 0, '<': 84, '<?xml version=': 0, '>': 84, 'https://': 6, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '110591', 'license': 'MIT', 'n_tokens_mistral': 7609, 'n_tokens_neox': 7247, 'n_words': 2995, 'path': 'asteroid_notebook_Walter_20180006.ipynb', 'repo_name': 'hwalter734/FINALPROJECT-EOML', 'size': '23664'}