Usamos un catálogo de galaxias identificadas en VVV en los tiles d010 y d0115 de Baravalle L.
Para saber donde estan ubicados los tiles usamos el mapa de VVV
En estos tiles encontraron 574 objetos con propiedades morfologicas, fotometricas y fotocromaticas propias de galaxias. 90 de los mismos han sido visualmente inspeccionados, y constituyen una muestra bona fide de galaxias en el VVV.
In [17]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display
from astropy.io import ascii
from astropy.table import Table, Column
%matplotlib inline
Leemos los datos de la tabla de Baravalle et al. 2017
El catálogo se encuentra en formato cds, pero hubo que modificar los datos perdidos.
Simplemente se reemplazaron los valores -- por espacios vacios.
In [66]:
cat = ascii.read('Table1.all.txt', format='cds')
Se crea una columna donde podamos marcar los objetos identificados visualmente.
In [67]:
visual = []
for row in cat:
if row['Id'].endswith('*'):
visual.append(1)
else: visual.append(0)
In [68]:
vis = Column(data=visual, name='Visual', dtype='Bool')
cat.add_column(vis)
El catalogo entonces quedaría con esta pinta:
In [71]:
cat
Out[71]:
Guardamos el catalogo en formato csv para conservar la columna con identificacion visual
In [72]:
cat.write('cat_lau.csv', format='csv', overwrite=True)
Veamos las distribuciones y dependencias de los datos.
In [73]:
pdcat = cat.to_pandas()
print pdcat.columns
In [74]:
numeric_cols = pdcat[[u'Z', u'Y', u'J', u'H', u'Ks', u'Z$_{2\prime\prime}$', u'Y$_{2\prime\prime}$',
u'J$_{2\prime\prime}$', u'H$_{2\prime\prime}$', u'Ks$_{2\prime\prime}$',
u'R$_{1/2}$', u'C', u'$\epsilon$', u'n']]
pd.scatter_matrix(numeric_cols, alpha=0.3, figsize=(12, 12), diagonal='kde')
plt.show()
Ahora que hemos chequeado las distribuciones de los valores de magnitudes y parametros extra podemos ver que a simple vista la muestra proviene de una misma poblacion de objetos.
In [103]:
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn import svm
from sklearn import neighbors
import sklearn.cross_validation as cv
from sklearn.model_selection import StratifiedKFold
from sklearn import metrics
In [76]:
def experiment(clf, x, y, nfolds=10):
skf = StratifiedKFold(n_splits=nfolds)
probabilities = np.array([])
predictions = np.array([])
y_testing = np.array([])
for train, test in skf.split(x, y):
x_train = x[train]
y_train = y[train]
clf.fit(x_train, y_train)
x_test = x[test]
y_test = y[test]
pr = clf.predict(x_test)
probs = clf.predict_proba(x_test)[:, 0]
probabilities = np.hstack([probabilities, probs])
predictions = np.hstack([predictions, pr])
y_testing = np.hstack([y_testing, y_test])
print metrics.classification_report(y_testing, predictions)
fpr, tpr, thresholds = metrics.roc_curve(y_testing, 1.-probabilities)
roc_auc = metrics.auc(fpr, tpr)
return {'fpr': fpr,
'tpr': tpr,
'thresh': thresholds,
'roc_auc': roc_auc,
'y_test': y_testing,
'predictions': predictions,
'probabilities': probabilities,
'confusion_matrix': metrics.confusion_matrix(y_testing, predictions),
}
In [98]:
X = pdcat[[u'J', u'H', u'Ks',
u'J$_{2\prime\prime}$', u'H$_{2\prime\prime}$', u'Ks$_{2\prime\prime}$',
u'R$_{1/2}$', u'C', u'$\epsilon$', u'n']]
X = X.as_matrix()
# Y = pdcat[u'Visual']
# Y = Y.as_matrix()
Y = np.asarray(pdcat['Visual'], dtype="int")
In [99]:
X = np.delete(X, [[378, 442]], axis=0)
Y = np.delete(Y, [[378, 442]], axis=0)
In [121]:
res_Dtree = experiment(DecisionTreeClassifier(max_leaf_nodes=5), X, Y, 10)
In [122]:
res_knn = experiment(neighbors.KNeighborsClassifier(n_neighbors=3, weights='uniform'), X, Y, 10)
In [123]:
res_svc = experiment(svm.SVC(kernel='linear', probability=True), X, Y, 10)
In [124]:
res_Rforest = experiment(RandomForestClassifier(n_estimators=100, min_samples_leaf=5), X, Y, 10)
In [125]:
res_C45 = experiment(DecisionTreeClassifier(criterion='entropy', max_leaf_nodes=5), X, Y, 10)
Ahora que probamos con varios clasificadores podemos construir una curva ROC
In [127]:
plt.figure(figsize=(6, 6))
#plt.figaspect(.8)
fpr = res_Dtree['fpr']
tpr = res_Dtree['tpr']
roc_auc = metrics.auc(fpr, tpr)
plt.plot(fpr, tpr, lw=1.5, color='blue', label='Dtree AUC={:06.4f}'.format(roc_auc))
fpr = res_C45['fpr']
tpr = res_C45['tpr']
roc_auc = metrics.auc(fpr, tpr)
plt.plot(fpr, tpr, lw=1.5, color='red', label='C45 AUC={:06.4f}'.format(roc_auc))
fpr = res_Rforest['fpr']
tpr = res_Rforest['tpr']
roc_auc = metrics.auc(fpr, tpr)
plt.plot(fpr, tpr, lw=1.5, color='green', label='Rforest AUC={:06.4f}'.format(roc_auc))
fpr = res_svc['fpr']
tpr = res_svc['tpr']
roc_auc = metrics.auc(fpr, tpr)
plt.plot(fpr, tpr, lw=1.5, color='black', label='SVM AUC={:06.4f}'.format(roc_auc))
fpr = res_knn['fpr']
tpr = res_knn['tpr']
roc_auc = metrics.auc(fpr, tpr)
plt.plot(fpr, tpr, lw=1.5, color='magenta', label='KNN AUC={:06.4f}'.format(roc_auc))
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.legend(loc='lower right')
plt.title('Receiver Operating Characteristic Curve')
Out[127]:
In [ ]: