<jupyter_start><jupyter_text>Exercises and supplementary notes - Week 4, Session AM1We'll be using the same tools that we used last week for this session.- [pandas](pandas.pydata.org) for data handling (our dataframe library)- [seaborn](seaborn.pydata.org) for _nice_ data visualization- [scikit-learn](scikit-learn.org) an extensive machine learning library.- [numpy](numpy.org) - a fundamental maths library best used by people with a strong maths background.  We won't explore it much today, but it does have some useful methods that we'll need.  It underlies all other mathematical and plotting tools that we use in Python.Shortly we'll also by trying out:- [statsmodel](statsmodel.org) - this is another library for doing statistical fitting.  It generates R-like reports.- We'll also being trying out two new parts of scikit-learn sklearn.cross_decomposition and decomposition for PLS and PCA We'll be using scikit-learn over the next few weeks, and it's well worth reading the documentation and high level descriptions.As before, the aim is to get familiar with code-sharing workflows - so we will be doing pair programming for the duration of the day! _You will probably want to take a moment to look at the documentation of the libraries above - especially pandas_The other useful resource is Stack Overflow - if you have a question that sounds like 'how do I do {x}' then someone will probably have answered it on SO. Questions are also tagged by library so if you have a particular pandas question you can do something like going to https://stackoverflow.com/questions/tagged/pandas (just replace the 'pandas' in the URL with whatever library you're trying to use.Generally answers on SO are probably a lot closer to getting you up and running than the documentation. Once you get used to the library then the documentation is generally a quicker reference. We will cover strategies for getting help in class.Topics that we'll be discussing in this session include:Robust Regression - http://scikit-learn.org/stable/modules/linear_model.htmlrobustness-regression-outliers-and-modeling-errors Git linksWe will be working through using GitHub and GitKraken to share code between pairs. We will go through all the workflow in detail in class but here are some useful links for reference:- GitKraken interface basics: https://support.gitkraken.com/start-here/interface- Staging and committing (save current state -> local history): https://support.gitkraken.com/working-with-commits/commits- Pushing and pulling (sync local history  GitHub history): https://support.gitkraken.com/working-with-repositories/pushing-and-pulling- Forking and pull requests (request to sync your GitHub history  someone else's history - requires a _review_):  - https://help.github.com/articles/about-forks/  - https://help.github.com/articles/creating-a-pull-request-from-a-fork/  Exercise:  Apply statsmodel to the synthetic drilling hole data.Statsmodel has an API with similarities to scikit-learn, but uses statstical language (particularly as used in financial and economic models) rather than the terminology that is more common in machine learning.  Statsmodel refers to endogeneous and exogeneous variables.  In many ways they reflect the differences in philosophy between how people with a statistics modelling background work, vs people with machine learning/computing backgrounds.  Scikit-learn has a focus on training and validation error curves and cross-validation to choose a model, whereas statsmodel provides metrics for hypothesis tests and goodness-of-fit.We'll briefly look at a typical report that statsmodel generates after fitting.<jupyter_code>import statsmodels.api as sm
import numpy as np
import matplotlib.pyplot as plt

from week4_help import create_composition_dataframe, create_observations, create_templates_matrix

templates = create_templates_matrix()
compositions = create_composition_dataframe(150)
observations = create_observations(compositions, templates)
plt.plot(observations[:, 3])
X = templates
y = observations[:, 3]
model = sm.OLS(y, X).fit()
model.summary()<jupyter_output><empty_output><jupyter_text>Let's compare this against truth:<jupyter_code>compositions.head()<jupyter_output><empty_output><jupyter_text>Exercise Look at the correlation matrixLets look at the correlation matrix for the features.  As you've seen, the feature variables are highly correlated with each other.  In real-world situations with more noise than this synthetic data this can cause linear least squares to fail.  Shortly we'll look at dimension reduction in the case of highly correlated variables, but lets look at the correlation visually.Seaborn's heatmap is useful for this, but it can be very slow.  Matplotlib's matshow is faster, try this if heatmap doesn't work well on your machine.<jupyter_code>import seaborn as sns
correlation_matrix = np.corrcoef(X)
sns.heatmap(correlation_matrix)
plt.matshow(correlation_matrix)<jupyter_output><empty_output><jupyter_text>Exercise:  Apply PCA to observations matrix.<jupyter_code>from sklearn.decomposition import PCA
pca = PCA(n_components=15)
pca.fit(observations)
print(pca.components_)<jupyter_output>[[ 0.04356901  0.04441032  0.04486882 ...  0.08874902  0.08618533
   0.09328089]
 [ 0.14093176  0.14811396  0.1402265  ...  0.06585451  0.05540338
  -0.06034275]
 [ 0.03005075  0.08117669  0.08143129 ... -0.06514712 -0.13922057
   0.05005523]
 ...
 [ 0.02451159 -0.04774677  0.05408418 ... -0.03501358 -0.06885206
   0.00726261]
 [-0.03701721  0.01998606 -0.12339401 ...  0.11767886  0.01833735
   0.08830087]
 [ 0.08597955 -0.05745995  0.01891607 ...  0.00737672 -0.09656674
   0.08652666]]<jupyter_text>Oops!  Can you see what's wrong here?  The components should have the same length as the orginal observations.  These components have only 50 elements.  This means that the observations array isn't organised properly for PCA fitting.  It thinks that we have ~400 observations, each with 50 features.  To fix this we should transpose the matrix.<jupyter_code>transposed_observations = observations.transpose()
pca.fit(transposed_observations)
print(pca.components_)<jupyter_output>[[-5.78987874e-04 -5.13463238e-04 -4.12026292e-04 ... -4.39840618e-04
  -2.21400840e-04  1.24469178e-04]
 [-3.20204687e-03  1.52245623e-03 -7.89177877e-04 ... -1.21231112e-03
  -3.30263693e-03 -3.67739212e-03]
 [ 1.78277824e-02 -5.87021932e-03  1.06311315e-04 ... -1.60123686e-03
   2.81856041e-03 -1.79621843e-04]
 ...
 [ 3.89525632e-03  3.96883291e-03 -1.40486631e-02 ...  3.41528003e-02
   2.97872963e-02 -2.87862988e-02]
 [ 8.82041236e-03 -2.23163296e-02 -9.23304956e-02 ... -3.76783244e-04
  -1.17443066e-01  6.27055760e-02]
 [-2.28986030e-02 -1.01964331e-01 -2.14187132e-02 ...  3.71879922e-04
  -1.02654439e-02  5.05762043e-02]]<jupyter_text>That's better!  We clearly have many more elements in these component arrays now.Let's look at the explained variance.  Warning:  Older documentation will tell you to do this. But Python is a living project and this now returns abstract quantities ("eigenvalues") that are related to the variances, but are not variances.<jupyter_code>print(pca.explained_variance_)<jupyter_output>[1.50295476e+05 1.95388466e+04 6.46458273e+02 2.67659009e+02
 6.20436882e+01 6.03735309e+01 5.86986570e+01 5.75539133e+01
 5.59864336e+01 5.49602111e+01 5.42284546e+01 5.33005057e+01
 5.23914273e+01 5.11941665e+01 5.06725681e+01]<jupyter_text>Nowadays we should use explained_variance_ratio, like this:<jupyter_code>print (pca.explained_variance_ratio_)<jupyter_output>[8.61734793e-01 1.12028016e-01 3.70653595e-03 1.53465085e-03
 3.55733959e-04 3.46157938e-04 3.36554874e-04 3.29991366e-04
 3.21004058e-04 3.15120104e-04 3.10924502e-04 3.05604010e-04
 3.00391714e-04 2.93527095e-04 2.90536456e-04]<jupyter_text>This says that the first component was able to account for about 86% of the variation, and the two following components account for about 1% and then less than 1% each.PCA has compressed 89% of the variation in ~400 features into just two transformed features!<jupyter_code>plt.plot(pca.components_[0])
plt.plot(templates[:,0])
plt.plot(pca.components_[1])
plt.plot(templates[:,1])
plt.plot(templates[:,2])
plt.plot(pca.components_[2])
plt.plot(templates[:, 3])
plt.plot(pca.components_[3])<jupyter_output><empty_output><jupyter_text>This looks better than I expected, for unobtainium.  But there's a good chance that PCA regression will struggle to predict unobtainium well, especially when there isn't very much present.We can plot explained variance like this:<jupyter_code>plt.plot(np.cumsum(pca.explained_variance_ratio_))
plt.xlabel('number of components')
plt.ylabel('cumulative explained variance')<jupyter_output><empty_output><jupyter_text>We seem to hit a sweet spot at about 2 components, and then explained variance only gradually improves after this point.  This is a good result, given that we know that we know these samples only contain four minerals, and all other variance is from noise.How do we look at the transformed data?  We use the transform method to find the "scores" associated with each observation.  These are new observations in the new feature space that PCA has found.  Each element in observations_in_pc_space has 15 elements in it, being the reduced feature space.<jupyter_code>observations_in_pc_space = pca.transform(transposed_observations)
plt.plot(observations_in_pc_space[0])
plt.plot(observations_in_pc_space[2])<jupyter_output><empty_output><jupyter_text>We'd like to be able to plot the scores against each other.  Let's plot PC1 against PC2 for all observations.<jupyter_code>plt.plot(observations_in_pc_space[:, 0],
        observations_in_pc_space[:,1],
        'o')<jupyter_output><empty_output><jupyter_text>Oh look!  There are clusters!  What about PC2 and PC3?<jupyter_code>plt.plot(observations_in_pc_space[:, 1],
        observations_in_pc_space[:,2],
        'o')<jupyter_output><empty_output><jupyter_text>Try to study the Octane dataset in the same manner!  You'll need to disregard the column with the octane ratings, as we're only interested in the spectra.<jupyter_code>import os
from os import path
import pandas as pd

current_folder = os.getcwd()

print ('Current working folder: ' + current_folder)

file_name = 'octane.xlsx'
sheet_name = 'Sheet1'

df = pd.ExcelFile(path.join(current_folder, file_name)).parse(sheet_name)
#df.set_index('Sample name', inplace=True)
octane = df.drop(columns=['Octane number','Sample name'])
octane.head()
import seaborn as sns
import numpy as np
correlation_matrix = np.corrcoef(octane.transpose())
sns.heatmap(correlation_matrix)
from sklearn.decomposition import PCA
pca = PCA(n_components=15)
pca.fit(octane)
print(pca.components_)
pca.fit(octane)
print(pca.explained_variance_)
print (pca.explained_variance_ratio_)
import matplotlib.pyplot as plt
plt.plot(pca.components_[0])
plt.plot(np.cumsum(pca.explained_variance_ratio_))
plt.xlabel('number of components')
plt.ylabel('cumulative explained variance')
observations_in_pc_space = pca.transform(octane)
plt.plot(observations_in_pc_space[0])
plt.plot(observations_in_pc_space[:, 0],
        observations_in_pc_space[:,1],
        'o')
plt.plot(observations_in_pc_space[:, 1],
        observations_in_pc_space[:,2],
        'o')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08766032728881026, 'fraction_numerical': 0.09774436090225563, 'mean_word_length': 5.151251360174102, 'pattern_counts': {'":': 0, '<': 49, '<?xml version=': 0, '>': 50, 'https://': 6, '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': '13', 'dir': 'jupyter-structured-clean-dedup', 'id': '616841', 'license': 'MIT', 'n_tokens_mistral': 4307, 'n_tokens_neox': 3597, 'n_words': 1156, 'path': 'notebooks/week4_am2.ipynb', 'repo_name': 'TimDaviesATCO/04-multivariate-analysis', 'size': '11305'}