<jupyter_start><jupyter_text>**Exercise set 5**==================>The goal of this exercise is to perform **principal component analysis**,>**clustering**, and **classification** on a data set with many variables.**Exercise 5.1**In this exercise, we will explore the "wine dataset" which is a common exampledataset used for classification. The dataset contains the results ofa chemical analysis of wines from the same region in Italy, using grapes grownby three different cultivators. In this first exercise, we will explore thisdataset using principal component analysis.**(a)** Begin by exploring the raw data. Here, you should choosethe method yourself. You can, for instance, look at histograms of thedifferent measured quantities, correlations between the quantities,or other plots of the raw data. It can also be useful to explorestatistical properties like averages and standard deviations. The Pythoncode in the following cell can be used to load the data set, and it will print out some summaries of the raw data which you may findhelpfull for your explorationAfter looking at the raw data, are there some of the variables that seem to be able to distinguish between the wines produced by the different cultivators?<jupyter_code>"""Load the wine data set and print some info."""                               
from sklearn.datasets import load_wine                                          
import numpy as np                                                              
import pandas as pd                                                             
                                                                                
                                                                                
data_set = load_wine()                                                          
# Print out some information about the data set:                                
print('Variables in the data set:')                                             
for i in data_set['feature_names']:                                             
    print(i)                                                                    
print('Classes in the data set (cultivators):')                                 
for i in data_set['target_names']:                                              
    print(i)                                                                    
# Convert the data set into a pandas DataFrame:                                 
data = pd.DataFrame(data_set['data'], columns=data_set['feature_names'])        
# Print a table with a summary for each variable:                               
for variable in data.columns:                                                   
    print('\nSummary for "{}":'.format(variable))                               
    print(data[[variable]].describe())                                          
# We can also use the class information:                                        
class_data = data_set['target']                                                 
class_names = dict(enumerate(data_set['target_names']))                         
variable = 'color_intensity'                                                    
for class_id, class_name in class_names.items():                                
    print('\nInformation about "{}" for "{}"'.format(variable, class_name))        
    idx = np.where(class_data == class_id)[0]                                   
    data_class = data.loc[idx, variable]                                        
    print(data_class.describe())
# your code here<jupyter_output><empty_output><jupyter_text>**Your answer to 5.1(a)**: *double click here* **(b)** Perform a PCA on the data set and plot the explained variance as a functionof the number of principal components. Do you need to scale your data before performingPCA in this case (why/why not)? How many principal components are needed to explain $95$ \%of the variance in the data? The following code cell can be used to run the PCA<jupyter_code>"""Load the wine data set and run PCA."""                                       
from sklearn.datasets import load_wine                                          
from sklearn.preprocessing import scale                                         
from sklearn.decomposition import PCA                                           
import numpy as np                                                              
import pandas as pd                                                             
                                                                                
                                                                                
data_set = load_wine()                                                          
data = pd.DataFrame(data_set['data'], columns=data_set['feature_names'])        
X = data                                                                        
# Uncomment the following line to scale your data:                              
#X = scale(data)                                                                
pca = PCA()                                                                     
scores = pca.fit_transform(X)                                                   
# Print out the percentage of variance explained by each component:             
print(pca.explained_variance_ratio_)<jupyter_output><empty_output><jupyter_text>**Your answer to 5.1(b)**: *double click here* **(c)** After selecting the number of principal components to use, rerun the PCA withthis number of components. This can be done by defining argument `n_components` to `PCA`, e.g.: `pca = PCA(n_components=13)`.Obtain the scores, and make a plot of the scores forprincipal component 1 (on the x-axis) and principal component 2 (on the y-axis).Do you see any grouping(s) ("clusters") in your scores plot?Here, you can choose to color the scores accordingto the class they belong to (i.e. by using the classdata available in the data set).<jupyter_code># Your code here<jupyter_output><empty_output><jupyter_text>**Your answer to 5.1(c):** *(Double click here)* **(d)** Explore the loadings for your PCA model by plotting theloadings for the variables (on principal component 1 andprincipal component 2). Do any of the variables seem to be correlated?<jupyter_code># Your code here<jupyter_output><empty_output><jupyter_text>**Your answer to 5.1(d):** *(Double click here)* **(e)** Save the scores you have obtained to a new file.We will use this information in the next partof the exercise, where we will try to find clusters in our data.Saving the scores can be done with `pandas` as follows:<jupyter_code>"""Load the wine data set, run PCA and save scores."""
from sklearn.datasets import load_wine
from sklearn.preprocessing import scale
from sklearn.decomposition import PCA
import numpy as np
import pandas as pd


data_set = load_wine()
data = pd.DataFrame(data_set['data'], columns=data_set['feature_names'])
X = data
# Uncomment the following line to scale your data:
#X = scale(data)
pca = PCA()
scores = pca.fit_transform(X)
# Create variable names for the principal components:
pc_name = ['PC{}'.format(i + 1) for i in range(pca.n_components_)]
# Create a DataFrame from the scores:
scores_data = pd.DataFrame(scores, columns=pc_name)
# Save the scores to a comma separated values-file:
scores_data.to_csv('scores.csv')<jupyter_output><empty_output><jupyter_text>**Exercise 5.2**We will continue exploring the "wine dataset". We will pretend that we do notknow that there are $3$ classes in the dataset, and we will investigatewhat the `KMeans` clustering method can tell us about it. For thisexercise, it is a good idea to read through all points below beforestarting, as you will find a link to a specific example you can use to answer most of the questions.**(a)** Explain the steps in the `KMeans` clustering algorithm.How can we use this algorithm without knowing how many clustersthere are in the data? **Your answer to 5.2(a):** *(Double click here)* **(b)**Run `KMeans` clustering on the wine dataset. Here, you will have toselect a set of numbers of clusters to look for. (Limit yourself toa maximum of 10 clusters) After running the clustering for your data, plot the sum of squared distances of samples to their closestcluster center, as a function of the number of clusters considered. Explain briefly how this plot can be used to identify the "correct"number of clusters. How many clusters would you say there are in thedataset, based on this plot alone?      To get you started, the cell below has some Python code that can be used to run theclustering and store the results (see also the [silhouette example](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.silhouette_samples.html))Note that the `cluster_km` object contains the following results as attributes: * `cluster_centers_`: Coordinates of cluster centers. * `labels_`: Labels of each point. * `inertia_`: Sum of squared distances of samples to their closest cluster center. * `n_iter_`: Number of iterations run.<jupyter_code>"""Load the wine data set and run KMeans."""                                    
from sklearn.datasets import load_wine                                          
from sklearn.preprocessing import scale                                         
from sklearn.cluster import KMeans                                              
import numpy as np                                                              
import pandas as pd                                                             
                                                                                
                                                                                
data_set = load_wine()                                                          
data = pd.DataFrame(data_set['data'], columns=data_set['feature_names'])        
X = scale(data)                                                                 
# Define a set of numbers of clusters to run KMeans for:                        
number_of_clusters = [2, 5]                                                     
# Set up variables for storing the results                                      
results = []  # Results for the clustering                                      
yfit = []  # Predicted clusters for data points in X                            
for i in number_of_clusters:                                                    
    cluster_km = KMeans(                                                        
        n_clusters=i,                                                           
        init='k-means++',                                                       
    )                                                                           
    y = cluster_km.fit_predict(X)                                               
    results.append(cluster_km)                                                  
    yfit.append(y)                                                              
# Print out some results:                                                       
print('Sum of squared distances of samples to their closest cluster center:')   
for i, result in zip(number_of_clusters, results):                              
    print('Clusters: {}: {}'.format(i, result.inertia_))
# Your code here<jupyter_output><empty_output><jupyter_text>**Your answer to 5.2(b):** *(Double click here)* **(c)** A general method that can be used to assess the clustering, is the silhouette method. This method calculates a silhouette value for each object which is a measure of how similar the object is to the cluster it belongs to (cohesion) compared toother clusters (separation). This is rather easy to calculate with `sklearn` as there is a method to do just so: [`silhouette_samples`](https://scikit-learn.org/stable/modules/generated/sklearn.metrics.silhouette_samples.html) from the module `sklearn.metrics`.Do the following: * For each clustering you have considered, i.e. for each number of clusters you tried,   calculate the silhouette values. * Plot the average silhouette value as a function of the number of clusters considered. * For each clustering, plot the silhouette values grouped into clusters. Say, if you,   for instance, considered $4$ clusters in one of your clusterings, plot the silhouette   values for each of these $4$ clusters. An example of how to do this is available on   the website of [`sklearn`](https://scikit-learn.org/stable/auto_examples/cluster/plot_kmeans_silhouette_analysis.html).   Using these results (average silhouette values) and the plots of silhouette values, whatis the best number of clusters to use? How does this compare with what we already know --that the samples come from $3$ different cultivators?<jupyter_code># Your code here<jupyter_output><empty_output><jupyter_text>**Your answer to 5.2(c):** *(Double click here)* **(d)** Rerun your analysis on the scores you stored in the last point of the PCA part.But use only the scores from principal components $1$ and $2$.Do the results from this analysis differ from the cluster analysis on the full data set?**Note:** As we only consider two of the principal components here, we have 2D-data. Thismeans that we can plot the clusters more easily. If you are curious, plot thescores for principal components $1$ and $2$ and color the points according to theclustering results you have obtained. Here, you can also show the centers of theclusters by using the `cluster_centers_` attribute of the `KMeans` object you haveused. This part of the exercise also shows that PCA can be used as an initial method to reduce the dimensionality of the original problem. We have here combined PCA and KMeans to solve a clustering problem.<jupyter_code># Your code here<jupyter_output><empty_output><jupyter_text>**Your answer to 5.2(d):** *(Double click here)* **Exercise 5.3: LDA Example**Both PCA and KMeans are examples of unsupervised methods -- we did not use the class information available to us tofind clusters in our data. [Linear discriminant analysis](https://en.wikipedia.org/wiki/Linear_discriminant_analysis)on the other hand, is a supervised method that uses the classinformation for *classification*. In your own words, how would youdescribe the difference between *classification* and*clustering*?LDA is similar to PCA, but rather than looking for latent variables that maximizethe covariance in our data, we rather look for latent variables that maximize the *class separation*.Below, you will find a small script that will run LDA on thewine data set. Run this script and observe the results.Note here the difference when we train the LDAmodel: `X_trans = lda.fit_transform(X, y)`.We are supplying the "y" values (i.e. the classes) which is what we expectfor a supervised method.For the curious student: Apply LDA to the 2D-example dataset from exercise 4,where we investigated classification by PCA. Does this classification differ fromthe simple rule we found there?<jupyter_code>"""Load the wine data set and run LDA."""                                       
from matplotlib import pyplot as plt                                            
from matplotlib.colors import ListedColormap                                    
from matplotlib.cm import tab10                                                 
from sklearn.datasets import load_wine                                          
from sklearn.preprocessing import scale                                         
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis            
import numpy as np                                                              
import pandas as pd                                                             
plt.style.use('seaborn-talk')                                                   
                                                                                
                                                                                
data_set = load_wine()                                                          
data = pd.DataFrame(data_set['data'], columns=data_set['feature_names'])        
X = scale(data)                                                                 
y = data_set['target']  # load the class information                            
# Run LDA:                                                                      
lda = LinearDiscriminantAnalysis()                                              
X_trans = lda.fit_transform(X, y)                                               
print('Number of classes:', len(lda.classes_))                                  
# Predict classes for our original points:                                      
y_hat = lda.predict(X)                                                          
                                                                                
# Plot the explained variance:                                                  
fig1, ax1 = plt.subplots()                                                      
comp = list(range(1, len(lda.explained_variance_ratio_) + 1))                   
ax1.bar(comp, lda.explained_variance_ratio_,                                    
        label='Variance explained by component')                                
ax1.plot(                                                                       
    [0] + comp,                                                                 
    [0] + list(np.cumsum(lda.explained_variance_ratio_)),                       
    color='black',                                                              
    marker='o',                                                                 
    label='Cumulative variance explained')                                      
ax1.set_xticks([0] + comp)                                                      
ax1.set(xlabel='LDA component', ylabel='Ratio of variance explained')           
ax1.legend()                                                                    
ax1.axhline(y=1, ls=':', color='black', alpha=0.8)                              
fig1.tight_layout()                                                             
                                                                                
# Plot the transformed X, this is similar to the scores found in PCA:           
fig2, ax2 = plt.subplots()                                                      
for i in np.unique(y_hat):                                                      
    ax2.scatter(X_trans[y_hat == i, 0], X_trans[y_hat == i, 1],                 
                color=tab10.colors[i], s=150)                                   
ax2.set(xlabel='LDA component 1', ylabel='LDA component 2')                     
# Plot the centers of the clusters found:                                       
for center in lda.transform(lda.means_):                                        
    ax2.scatter(center[0], center[1], s=250, color='black',                     
                marker='X', edgecolor='white')                                  
                                                                                
# Now, in order to plot the regions, we would like to have 2D data.             
# The classification we have right now, expects 13 variables to                 
# classify samples. We therefore run a second LDA on the LDA we                 
# already have performed:                                                       
lda2 = LinearDiscriminantAnalysis()                                             
X_trans2 = lda2.fit_transform(X_trans, y_hat)                                   
y_hat2 = lda2.predict(X_trans)                                                  
# Show the regions:                                                             
fig3, ax3 = plt.subplots()                                                      
X_set = X_trans
X1, X2 = np.meshgrid(                                                           
    np.linspace(X_trans[:, 0].min() - 1, X_trans[:, 0].max() + 1, 500),         
    np.linspace(X_trans[:, 1].min() - 1, X_trans[:, 1].max() + 1, 500)          
)                                                                               
Z = lda2.predict(np.array([X1.ravel(), X2.ravel()]).T).reshape(X1.shape)        
plt.contourf(X1, X2, Z, alpha = 0.3,                                            
             cmap=ListedColormap(tab10.colors[:3]))                             
# Add the original samples:                                                     
for i in np.unique(y_hat2):                                                     
    ax3.scatter(X_trans2[y_hat2 == i, 0], X_trans2[y_hat2 == i, 1],             
                color=tab10.colors[i], s=150)                                   
ax3.set(xlabel='LDA component 1', ylabel='LDA component 2')                     
# Plot the centers of the clusters found:                                       
for center in lda2.transform(lda2.means_):                                      
    ax3.scatter(center[0], center[1], s=250, color='black',                     
                marker='X', edgecolor='white')                                  
fig3.tight_layout()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.0644674835061263, 'fraction_numerical': 0.006974552309142319, 'mean_word_length': 1.3620881567230632, 'pattern_counts': {'":': 1, '<': 37, '<?xml version=': 0, '>': 39, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 191, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '133498', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 4599, 'n_tokens_neox': 4312, 'n_words': 1825, 'path': 'exercises_2020/05_Exercise_Set_5.ipynb', 'repo_name': 'sroet/chemometrics', 'size': '21220'}