cas-pml/NSL/aufgaben/notebooks/solutions.ipynb

Solution Exercise 1: k-Means Clustering Unsupervised Learning CAS Practical Machine Learning, SS-HS23 Prof. Dr. Matthias Dehmer

In [1]:

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

Load and Prepare Dataset

Load data from tab-separated file using Pandas and show first rows.

In [2]:

df = pd.read_csv('data/data_delivery_fleet.tsv', sep='\t', header=0)
df.head()

Out [2]:

	Driver_ID	Distance_Feature	Speeding_Feature
0	3423311935	71.24	28.0
1	3423313212	52.53	25.0
2	3423313724	64.54	27.0
3	3423311373	55.69	22.0
4	3423310999	54.58	25.0

Column 1: Driver_ID is the identification of the delivery driver.
Column 2: Distance_Feature is the average number of miles driven per day.
Column 3: Speeding_Feature is the % of time the driver exceeded the speed limit by more >5 mph.

Create NumPy array data from the 2nd (Distance_Feature) and the 3rd (Speeding_Feature) column.

In [8]:

data = df[['Distance_Feature', 'Speeding_Feature']].as_matrix()
print(data.shape)
print(data.__class__)
data

Out [8]:

(4000, 2)
<class 'numpy.ndarray'>

array([[  71.24,   28.  ],
       [  52.53,   25.  ],
       [  64.54,   27.  ],
       ..., 
       [ 170.91,   12.  ],
       [ 176.14,    5.  ],
       [ 168.03,    9.  ]])

In [4]:

plt.scatter(df.Distance_Feature, df.Speeding_Feature, marker='.')
plt.xlabel('distance [mi]')
plt.ylabel('speeding fraction [%]')
plt.grid()
plt.show()

Cluster the Data

Initially, pick some cluster count k that you see fit.

In [9]:

from sklearn.cluster import KMeans

k_means = KMeans(n_clusters=4).fit(data)

Visualize the clusters.

In [10]:

def visualize_k_means(data, k_means):
    plt.figure()
    k = len(k_means.cluster_centers_)
    for i in range(k):
        colors = ['r', 'g', 'b', 'y', 'c', 'm']
        points_in_cluster = np.array([data[j] for j in range(len(data)) if k_means.labels_[j] == i])
        plt.scatter(points_in_cluster[:, 0], points_in_cluster[:, 1], \
                    marker='.', color=colors[i % len(colors)])
    plt.scatter(k_means.cluster_centers_[:,0], k_means.cluster_centers_[:, 1], marker='*', color='k')
    plt.title('k-Means Clustering with k={} clusters'.format(k))
    plt.show()
    
visualize_k_means(data, k_means)

Plot cluster assignments for k=2 ... 8

In [11]:

for k in range(1,9):
    k_means = KMeans(n_clusters=k).fit(data)
    visualize_k_means(data, k_means)

Determine the optimal number of clusters

What is the best number of clusters?

Within-Cluster Sum of Squares:

In [8]:

Ks = range(1, 10)
inertia = []
for k in Ks:
    k_means = KMeans(n_clusters=k, random_state=0).fit(data)
    inertia.append(k_means.inertia_)
plt.bar(Ks, inertia)
plt.title('Within-Cluster Sum of Squares')
plt.xlabel('Cluster Count k')
plt.ylabel('Error')
plt.grid()
plt.show()

Percentage of Variance Explained

In [9]:

from scipy.spatial.distance import cdist, pdist

ks = range(1,10)  # Range of cluster count k 

# Fit the kmeans model for each n_clusters = k
k_means_var = [KMeans(n_clusters=k).fit(data) for k in ks]

# Pull out the cluster centers for each model
centroids = [X.cluster_centers_ for X in k_means_var]

# Calculate the Euclidean distance from 
# each point to each cluster center
k_euclid = [cdist(data, cent, 'euclidean') for cent in centroids]
dist = [np.min(ke,axis=1) for ke in k_euclid]

# Total within-cluster sum of squares
wcss = [sum(d**2) for d in dist]

# The total sum of squares
tss = sum(pdist(data)**2)/data.shape[0]

# The between-cluster sum of squares
bss = tss - wcss
# This follows from the fact that tss = within-cluster sum of squares + between-cluster sum of squares.

# elbow curve
plt.plot(ks, bss/tss*100, 'b*-')
plt.ylim((0,100))
plt.grid(True)
plt.xlabel('Cluster Count k')
plt.ylabel('Percentage of variance explained')
plt.title('Variance Explained vs. k')
plt.show()

Shilouette Score:

In [10]:

from sklearn.metrics import silhouette_score

s = []
ks = range(2,10)
for k in ks:
    k_means = KMeans(n_clusters=k)
    k_means.fit(data)
    s.append(silhouette_score(data, k_means.labels_, metric='euclidean'))

plt.plot(ks, s)
plt.ylabel("Silhouette Score")
plt.xlabel("Cluster Count k")
plt.title("Silhouette for K-means")
plt.show()

Hartigan Rule:

In [12]:

def hartigan_K(data, threshold=12):
    """
    'list_of_tuples' is a list containing the points you want to cluster
    'threshold' optimizes goodness of fit values
    Hartigan recommends setting threshold to 10, but Chiang & Mirkin confirm up to 12
    returns integer, "correct" number of clusters
    """ 
    inertia_list = np.zeros(len(data) + 1) # initializes for maximum possible clusters
    num = 0                                          # counter
    H_Rule = threshold + 1                           # simply initializes above threshold to meet 'while' condition
    
    # NOTE: 'inertia' is equivalent to the sum of within-cluster distances to centroids
    
    while num < len(data) and H_Rule > threshold:
        kmn = KMeans(n_clusters = num + 1)
        kmn.fit(data)
        inertia_list[num + 1] = kmn.inertia_
        if num > 0:
            H_Rule = ((float(inertia_list[num]) / inertia_list[num + 1]) - 1) * (len(data) - num - 1)
            print('k={}, H={}'.format(num, H_Rule))
        num += 1
    
    if H_Rule > threshold:
        num += 1
    # NOTE: if while-loop reaches the number of K-Means clusters equal to the length of list_of_tuples
    # without hitting the threshold, then function returns trivial solution that there are N clusters
    # (where N is the number of points under observation)
    
    return num - 1 

In [15]:

k_h = hartigan_K(data, threshold=12)

k=1, H=33006.98608114748
k=2, H=1303.7793604071974
k=3, H=1516.1728593448447
k=4, H=1381.735661629919
k=5, H=1733.4844033305608
k=6, H=662.1126967486157
k=7, H=350.75237466138384
k=8, H=617.1329342659591
k=9, H=421.1173631740277
k=10, H=406.5804351771144
k=11, H=416.8955261064009
k=12, H=362.84126685539223
k=13, H=268.0986265446576
k=14, H=376.79250950835444
k=15, H=248.9447268836972
k=16, H=226.38879767516693
k=17, H=185.65530642910906
k=18, H=191.8693177888545
k=19, H=145.59881362029617
k=20, H=210.95991507388356
k=21, H=161.5512299419221
k=22, H=195.21709743460372
k=23, H=156.42882039212276
k=24, H=125.93954025322607
k=25, H=167.8606728704937
k=26, H=119.8879408188638
k=27, H=223.48717969820447
k=28, H=65.44257204515195
k=29, H=162.37657240138034
k=30, H=170.8192731833724
k=31, H=65.35696340152566
k=32, H=195.77143090288106
k=33, H=29.43858884956867
k=34, H=139.71302113473672
k=35, H=98.40681862560261
k=36, H=111.88434490461884
k=37, H=92.54939400233751
k=38, H=155.5924843864272
k=39, H=78.23527913880022
k=40, H=95.09664615266315
k=41, H=62.88245645849113
k=42, H=105.62113726990547
k=43, H=119.27139183204469
k=44, H=100.56588807097297
k=45, H=102.48092233385731
k=46, H=78.36494108635938
k=47, H=58.310525120606286
k=48, H=115.59511991518754
k=49, H=30.153525751253508
k=50, H=93.8432475985359
k=51, H=60.0538323805475
k=52, H=77.017005051429
k=53, H=57.70526530638853
k=54, H=97.30799415088985
k=55, H=146.5339327759247
k=56, H=20.280811743746234
k=57, H=35.055323990746615
k=58, H=111.63760008594201
k=59, H=73.46186425913484
k=60, H=37.70109121330057
k=61, H=71.58763162108372
k=62, H=15.50552866656769
k=63, H=94.09957511683425
k=64, H=84.99243443931567
k=65, H=35.62311902363756
k=66, H=120.78559143956336
k=67, H=19.197115200930405
k=68, H=78.51716023633182
k=69, H=45.70217913713537
k=70, H=35.61189066066217
k=71, H=87.27002119997172
k=72, H=35.80918849022347
k=73, H=71.90047454836504
k=74, H=67.647449779252
k=75, H=23.31903490098106
k=76, H=78.11795047032751
k=77, H=-11.542781367279975

In [14]:

k_h

Out [14]:

53

In [14]:

k_means = KMeans(n_clusters=k_h, random_state=0).fit(data)
visualize_k_means(data, k_means)

In [ ]: