<jupyter_start><jupyter_text>Introduction The NEO dataset from NASA's JPL, located - https://ssd.jpl.nasa.gov/sbdb_query.cgi, will be evaluated with unsupervised learning. The main objective of this analysis is to use clustering algorithms to find a clustering classification method that produces useful clustering information. One use, in particular, that would be valuable is to utilise clustering classification of for prospecting NEO satellites in space mining operations. The sheer number of potential minable bodies in our solar system is enormous, and we don't even have a significant part of the solar bodies cataloged. The fact that they are so far away and spread out makes prospecting until you find the right asteroid or comet to mine unfeasible. A clustering model can take a sample set of prospecting data and use it to classify clusters and tell the miners where to go prospect and mine without ever having been there before. This project has the potential to save trillions of dollars, or more, that would have otherwise have been spent on blind prospecting. The potential market value of the space mining industry is preposteriously large; just one asteroid that we know of 16-Pshche is estimated to be worth around $700 quadrillion. There are many types of NEO satalites like comets that are composed up of water, CO2, and other compounds and then theres several types of asteroids which can have varrious concentrations of valuable metals. It would be helpful to be able to know what all of these satellites are composed of for mining ooperations but prospecting them all is not fesable and that's what this project aims to correct; at least, virtually.<jupyter_code>import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA
from scipy.spatial.distance import pdist, squareform
from sklearn.cluster import DBSCAN
from sklearn.cluster import KMeans
from sklearn import datasets, metrics
from sklearn.metrics import davies_bouldin_score

%matplotlib inline<jupyter_output><empty_output><jupyter_text>Data Cleaning and Preparation<jupyter_code># Import data into DataFrame
df = pd.read_csv("NEO_Data.csv")
# Review dataset features
df.info()<jupyter_output><class 'pandas.core.frame.DataFrame'>
RangeIndex: 25391 entries, 0 to 25390
Data columns (total 75 columns):
 #   Column          Non-Null Count  Dtype  
---  ------          --------------  -----  
 0   id              25391 non-null  object 
 1   spkid           25391 non-null  int64  
 2   full_name       25391 non-null  object 
 3   pdes            25391 non-null  object 
 4   name            345 non-null    object 
 5   prefix          186 non-null    object 
 6   neo             25391 non-null  object 
 7   pha             25202 non-null  object 
 8   H               25196 non-null  float64
 9   G               13 non-null     float64
 10  M1              129 non-null    float64
 11  M2              49 non-null     float64
 12  K1              129 non-null    float64
 13  K2              49 non-null     float64
 14  PC              31 non-null     float64
 15  diameter        1267 non-null   float64
 16  extent          10 non-null     object 
 17  albedo          1208 non-null  [...]<jupyter_text>The features listed above in our data set are abbreviations for astronomical terms in relation to properties of NEO satallite observations and measurements. The definition and meaning of these properties and measurements will not be looked into initially for this proof of concept model analysis. However, model performance and cluster evaluation could benifit from an in-depth study and application of these variables to the training of clustering models. It is useful to be able to look at the variables from a black box perspective and produce a sucessful model. To look into the features full descriptions, the information is available at the dataset's source - https://ssd.jpl.nasa.gov/sbdb_query.cgi<jupyter_code># Inspect feature content
pd.set_option("max_columns", 80)
df.head(10)<jupyter_output><empty_output><jupyter_text>Below I'm dropping several features for various reasons. I drop many beacause they are sparsly populated and contaon mostly NaN data. Then I drop the catagorical data for a few reasons. One reason is because I'm trying to train a classifier and do not want to include bias in my results. Working with the catagorical and numerical data presents more challenges and limits my modeling options. Lastly, I'm trying to keep this initial investigation simple and basic to serve as a proof of concept.<jupyter_code># Drop features not useful for analysis
df = df.drop(
    [
        "class",
        "pha",
        "id",
        "spkid",
        "full_name",
        "neo",
        "pdes",
        "name",
        "prefix",
        "G",
        "M1",
        "M2",
        "K1",
        "K2",
        "PC",
        "extent",
        "GM",
        "IR",
        "H_sigma",
        "diameter_sigma",
        "orbit_id",
        "epoch",
        "epoch_mjd",
        "epoch_cal",
        "equinox",
        "tp",
        "tp_cal",
        "diameter",
        "albedo",
        "rot_per",
        "BV",
        "UB",
        "spec_B",
        "spec_T",
        "n_del_obs_used",
        "n_dop_obs_used",
        "two_body",
        "A1",
        "A2",
        "A3",
        "DT",
        "producer",
        "first_obs",
        "last_obs",
        "n_obs_used",
    ],
    axis=1,
)
# Review Dataframe again
df.info()
df.head(5)
df.isna().mean()
# Drop small number of nulls
df = df.dropna()
df.shape<jupyter_output><empty_output><jupyter_text>Removing rows with null values only eliminated 420 observations out of 25,391, or just 1.65% of the data. For this initial investigation, this is considered an acceptable loss of data and should not effect the results in a measurable way. In future development, methods of handling the null values will be implimented to maximize the integrity of the data.<jupyter_code># Check out the overall statistics
df.describe()
# Use for rapid model development
# df = df.sample(10000, random_state=42)
# df.shape<jupyter_output><empty_output><jupyter_text>Data Standardization<jupyter_code># Scale the data
scaler = StandardScaler()
scaled = scaler.fit_transform(df)
scaled = pd.DataFrame(scaled, columns=df.columns, index=df.index)
scaled.head()<jupyter_output><empty_output><jupyter_text>Dimension Reduction<jupyter_code># Dimension Reduction - Principal Component Analysis
pca = PCA(n_components=2)
pca_components = pca.fit_transform(scaled)
plt.figure(figsize=(10, 5))
plt.scatter(pca_components[:, 0], pca_components[:, 1])
plt.xticks([])
plt.yticks([])
plt.title("2D Principal Component Representation ")
plt.axis("off")
plt.show()<jupyter_output><empty_output><jupyter_text>It's difficult to find much cluster meaning in this 2D representation. PCA is not the best coice for visualization but it will still be very helpful with the clustering analysis. A furthur investigation recommendation is to try more visualization with t-SNE or UMAP. Even so, we can see there is some interpretable structure in the data points. This provides an indication that even by reducing all of the features down to two dimensions PCA is retaining enough data to identify clusters.    Next, I evaluate the performance of selecting different numbers of principal components so I can utilise dimension reduction most effectivly for the upcoming clustering.<jupyter_code>pca = PCA()
pca.fit(scaled)
# Variance explained by each component
pca.explained_variance_
# Variance explained by each component (as a percentage)
pca.explained_variance_ratio_
# Cumulative explained variance
np.cumsum(pca.explained_variance_ratio_)
# Scree plot
plt.plot(pca.explained_variance_ratio_, "-o")
plt.title("Scree Plot")
plt.xlabel("Number of Components")
plt.ylabel("Explained Varience Ratio")
plt.show()
plt.plot(np.cumsum(pca.explained_variance_ratio_), "-o")
plt.title("Scree Plot")
plt.xlabel("Number of Components")
plt.ylabel("Sum of Explained Variance Raitos")
plt.show()<jupyter_output><empty_output><jupyter_text>It looks like using first 10 principal components will account for the majority of the variance in the data.<jupyter_code>pca = PCA(n_components=10)
components = pca.fit_transform(scaled)

components = pd.DataFrame(
    components,
    columns=["pc1", "pc2", "pc3", "pc4", "pc5", "pc6", "pc7", "pc8", "pc9", "pc10"],
    index=scaled.index,
)
components.head()<jupyter_output><empty_output><jupyter_text>DBSCAN Clustering<jupyter_code># Perform initial clustering
clusters = DBSCAN(eps=0.8, min_samples=19, metric="euclidean", n_jobs=-1)
clusters.fit(components)
# Assign labels to the dataframe
df["label"] = clusters.labels_
df["label"] = "Cluster: " + df["label"].astype(str)
dbscan_clusters = df[df["label"] != "Cluster: -1"]
dbscan_clusters["label"].value_counts()
# Silhouette Score
labels = clusters.fit_predict(components)
metrics.silhouette_score(components, labels, metric="euclidean")
# Davies Bouldin Score
davies_bouldin_score(components, labels)<jupyter_output><empty_output><jupyter_text>Since we do not have a ground truth to compare clusters to at the moment, I will use the Silhouette Score and the Davies Bouldin Score for evaluating the clustering. The parameters for `eps` and `min_samples` will be optimised below.  DBSCAN Cluster Evaluation<jupyter_code>epsilon = 0.8
max_samples = 25
silhouettes = []
for samples in range(1, max_samples + 1):
    labels = DBSCAN(
        eps=epsilon, min_samples=samples, metric="euclidean", n_jobs=-1
    ).fit_predict(components)
    silhouette = metrics.silhouette_score(components, labels, metric="euclidean")
    silhouettes.append(silhouette)
plt.title("Silhouette Scores, eps = 0.8")
plt.xlabel("min_samples")
plt.ylabel("eps")
plt.plot(silhouettes, "-o")
print(silhouettes)<jupyter_output>[-0.2781448237914977, -0.503502760364474, -0.5354585896010304, -0.542694062423307, -0.5228178747165709, -0.49375973754971236, -0.48818675123279276, -0.4779169412693037, -0.44915041170440184, -0.43206379533320455, -0.3960304122276044, -0.39639921569018666, -0.3142908903483431, -0.29838875673624354, -0.37901917393455165, -0.4048345578615378, -0.31133006202204466, -0.2543621661615011, -0.2570466429825432, -0.2807080054850371, -0.30196882882858134, -0.3362280587537896, -0.30150818142391755, -0.3078552816544328, -0.30616302010423185]<jupyter_text>Test records of high performing models* eps = 0.5, min_samp = 5, sil = -0.4452* eps = 0.5, min_samp = 6, sil = -0.3547* eps = 0.5, min_samp = 7, sil = -0.2983* eps = 0.6, min_samp = 7, sil = -0.4447* eps = 0.6, min_samp = 8, sil = -0.3954* eps = 0.6, min_samp = 9, sil = -0.3734* eps = 0.6, min_samp = 10, sil = -0.3332* eps = 0.6, min_samp = 11, sil = -0.2269* eps = 0.7, min_samp = 15, sil = -0.3660* eps = 0.7, min_samp = 22, sil = -0.2893* eps = 0.7, min_samp = 23, sil = -0.1883* eps = 0.8, min_step = 14, sil = -0.2984* eps = 0.8, min_step = 18, sil = -0.2544* eps = 0.8, min_step = 19, sil = -0.2570* eps = 0.8, min_step = 23, sil = -0.3015* eps = 1, min_samp = 25, sil = 0.0618* eps = 1.1, min_samp = 22, sil = 0.1229* eps = 1.2, min_samp = 19, sil = 0.1716* eps = 1.3, min_samp = 12, sil = 0.2536* eps = 1.4, min_samp = 10, sil = 0.3241* eps = 1.5, min_samp = 9, sil = 0.4059* eps = 1.6, min_samp = 7, sil = 0.4059* eps = 1.7, min_samp = 6, sil = 0.5400* eps = 1.8, min_samp = 6, sil = 0.5888* eps = 1.9, min_samp = 6, sil = 0.6336* eps = 2.0, min_samp = 6, sil = 0.6653* eps = 2.1, min_samp = 7, sil = 0.6952* eps = 2.2, min_samp = 4, sil = 0.7424* eps = 2.3, min_samp = 4, sil = 0.7627* eps = 2.4, min_samp = 4, sil = 0.7960* eps = 2.5, min_samp = 4, sil = 0.8074* eps = 3.0, min_samp = 3, sil = 0.8074* eps = 3.5, min_samp = 4, sil = 0.8074 Now that the parameters have been optimised, as far as the silhouette score is concerned, it's time to take a look at the clusters produced by this set of parameters.<jupyter_code>params = [
    [0.8, 14],
    [0.8, 18],
    [0.8, 19],
    [0.8, 20],
    [0.8, 21],
    [0.8, 23],
]


""" Code for exploring clustering results using optimal parameters
    [0.5, 5],
    [0.5, 6],
    [0.5, 7],
    [0.6, 7],
    [0.6, 8],
    [0.6, 9],
    [0.6, 10],
    [0.6, 11],
    [0.7, 15],
    [0.7, 16],
    [0.7, 17],
    [0.7, 18],
    [0.7, 19],
    [0.7, 20],
    [0.7, 21],
    [0.7, 22],
    [0.7, 23],
    [0.8, 24],
    [0.8, 25],
    [0.8, 26],
    [0.8, 27],
    [0.8, 28],
    [0.8, 29],
    [0.8, 30],
    [1, 25],
    [1.1, 22],
    [1.2, 19],
    [1.3, 12],
    [1.4, 10],
    [1.5, 9],
    [1.6, 7],
    [1.7, 6],
    [1.8, 6],
    [1.9, 6],
    [2, 6],
    [2.1, 7],
    [2.2, 4],
    [2.3, 4],
    [2.4, 4],
    [2.5, 4],
"""
for param in params:
    clusters = DBSCAN(eps=param[0], min_samples=param[1], metric="euclidean", n_jobs=-1)
    clusters.fit(components)
    df["label"] = clusters.labels_
    df["label"] = "Cluster: " + df["label"].astype(str)
    dbscan_clusters = df[df["label"] != "Cluster: -1"]

    print("\n")
    print("Eps: " + str(param[0]) + ", min_samples: " + str(param[1]))
    print(dbscan_clusters["label"].value_counts())<jupyter_output>Eps: 0.8, min_samples: 14
Cluster: 1    4046
Cluster: 0      27
Cluster: 7      22
Cluster: 6      21
Cluster: 5      20
Cluster: 3      13
Cluster: 2      13
Cluster: 4       9
Cluster: 8       6
Name: label, dtype: int64


Eps: 0.8, min_samples: 18
Cluster: 0    2539
Cluster: 1     147
Cluster: 2      48
Cluster: 4      19
Cluster: 3      18
Cluster: 5      18
Name: label, dtype: int64


Eps: 0.8, min_samples: 19
Cluster: 0    2176
Cluster: 1     147
Cluster: 2     133
Cluster: 3      39
Cluster: 4      16
Cluster: 5      13
Name: label, dtype: int64


Eps: 0.8, min_samples: 20
Cluster: 0    1982
Cluster: 1     139
Cluster: 2     119
Cluster: 3      28
Cluster: 5      23
Cluster: 4      11
Name: label, dtype: int64


Eps: 0.8, min_samples: 21
Cluster: 0    1762
Cluster: 1     133
Cluster: 2      92
Cluster: 4      29
Cluster: 6      16
Cluster: 7      14
Cluster: 5      14
Cluster: 3      12
Name: label, dtype: int64


Eps: 0.8, min_samples: 23
Cluster: 0    1203
Cluster: 2     118[...]<jupyter_text>In nearly every case, the models with positive silhouette scores could only identify one cluster. Therefore, the models with the best negative silhouette scores should wil evaluated to select the best performing model. After reviewing the silhouette scores and clustering groups of several variations of combinations of eps and min_sample combinations, I find that using `eps = 0.8` and `min_sample = 19` to be one of the best classification models to be produced, resulting with `silhoette_score = -0.2570`, which is not terrible. It's a model with one of the closest scores to 0. It classifies 2524 satellites into 6 groups. True, 10% percent population classification seems poor, it can be improved, however it should still be sufficient to be useful for prospecting and mining, even at this stage of development.  K-means Clustering<jupyter_code># Fit kmeans
k = 6
clusters = KMeans(n_clusters=k, verbose=0, random_state=42)
clusters.fit(components)
df["label"] = clusters.labels_
df["label"] = "Cluster: " + df["label"].astype(str)

kmeans_clusters = df[df["label"] != "Cluster: -1"]
kmeans_clusters["label"].value_counts()<jupyter_output><empty_output><jupyter_text>K-means Clustering Evaluation<jupyter_code># Silhouette Score
labels = clusters.predict(components)
metrics.silhouette_score(components, labels, metric="euclidean")
# Davies Bouldin Score
davies_bouldin_score(components, labels)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10770214125918824, 'fraction_numerical': 0.08098434004474273, 'mean_word_length': 3.4486778504407165, 'pattern_counts': {'":': 0, '<': 55, '<?xml version=': 0, '>': 55, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, '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': '14', 'dir': 'jupyter-structured-clean-dedup', 'id': '663723', 'license': 'MIT', 'n_tokens_mistral': 6088, 'n_tokens_neox': 5205, 'n_words': 1795, 'path': 'NEO_Cluster_Analysis.ipynb', 'repo_name': 'SLoganQ/Unsupervised_Learning_NEO_Clustering_Analysis', 'size': '15645'}