NIR spectroscopy - Cookies

Require Python 3.7 Scikit-learn version ...

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from urllib.request import urlopen

from sklearn.model_selection import train_test_split
from sklearn import neighbors, tree

url="https://perso.univ-rennes1.fr/valerie.monbet/doc/cours/Biscuits.csv"
biscuits=np.loadtxt(urlopen(url),skiprows=1,delimiter=";")

Let us look at the spectra.

In [2]:
# Extraction de la colonne fat
fat=biscuits[:,0]
# Extraction des variables explicatives
X=biscuits[:,1:]
# Trace d'un spectre
plt.figure(1)
plt.plot(X[1,:])
plt.title("Un exemple de spectre")
# TIn the following plot, color is varying according to the fat percent
fat=fat/max(fat)
colors= plt.cm.inferno(fat)
plt.figure(2)
for i in range(len(fat)):
   plt.plot(X[i,:],color=colors[i])
plt.title("Spectres NIR")
plt.ylabel("Absorbances")
plt.show()

Nearest neighbors method

For each value of the nearest neighbor numbers, Monte Carlo cross validation is run with aout 25% of the samples in the validation set. At each iteration de Mean Square Error is stored.

First block with uniform weights.

In [3]:
from sklearn import neighbors

nb_neighbors = [1,5,10,20]


B = 100 # number of samples for the cross validation
mse_knn= np.zeros((B,len(nb_neighbors)))
for b in range(0,B):
    X_train, X_test, y_train, y_test = train_test_split(X,fat,test_size=9)
    mk=np.mean(X_train,axis=0)
    sk=np.maximum(np.std(X_train,axis=0),10*np.finfo(float).eps)
    X_train, X_test = np.add(X_train,-mk), np.add(X_test,-mk)
    X_train, X_test = np.multiply(X_train,1/sk),np.multiply(X_test,1/sk)
    for i,k in enumerate(nb_neighbors):
        knn = neighbors.KNeighborsRegressor(k, weights='uniform')
        knn.fit(X_train,y_train)
        y_pred = knn.predict(X_test)
        mse_knn[b,i] = np.mean((y_pred-y_test)**2)
    
plt.boxplot(np.sqrt(mse_knn),labels=nb_neighbors)
plt.ylabel('Root mean square error')
plt.xlabel('Nb of neighbors')

print("Best RMSE of knn method: ",np.round(np.sqrt(np.min(np.mean(mse_knn,axis=0))),2))

Best RMSE of knn method:  0.06

First block with inverse (euclidean) distance weights

In [4]:
mse_knn= np.zeros((B,len(nb_neighbors)))
for b in range(0,B):
    X_train, X_test, y_train, y_test = train_test_split(X,fat,test_size=9)
    mk=np.mean(X_train,axis=0)
    sk=np.maximum(np.std(X_train,axis=0),10*np.finfo(float).eps)
    X_train, X_test = np.add(X_train,-mk), np.add(X_test,-mk)
    X_train, X_test = np.multiply(X_train,1/sk),np.multiply(X_test,1/sk)
    for i,k in enumerate(nb_neighbors):
        knn = neighbors.KNeighborsRegressor(k, weights='distance')
        knn.fit(X_train,y_train)
        y_pred = knn.predict(X_test)
        mse_knn[b,i] = np.mean((y_pred-y_test)**2)
    
plt.boxplot(np.sqrt(mse_knn),labels=nb_neighbors)
plt.ylabel('Root mean square error')
plt.xlabel('Nb of neighbors')

print("Best RMSE of knn method: ",np.round(np.sqrt(np.min(np.mean(mse_knn,axis=0))),2))

Best RMSE of knn method:  0.05

In the examples above, the cross validation procedure is detailled. But, scikit-learn provide some tools to perform it.

Shuffle split generate the B training and validation sets.

In [5]:
from sklearn.model_selection import GridSearchCV, cross_val_score
from sklearn.model_selection import ShuffleSplit

# Generation of the B training/validation sets
cv = ShuffleSplit(n_splits=B, 
                  test_size=0.25, 
                  random_state=0)
# Initialization of the regressor
knn = neighbors.KNeighborsRegressor()
param_grid = [{'n_neighbors': [1,5,10,20], 
               'weights': ['uniform','distance'],
               'metric': ['euclidean','minkowski']}]
# Optimization of the hyperparameters
search = GridSearchCV(knn, 
                      param_grid, 
                      scoring="neg_mean_squared_error",
                      cv=cv,
                      n_jobs=-1)
search.fit(X,fat)
Out[5]:
GridSearchCV(cv=ShuffleSplit(n_splits=100, random_state=0, test_size=0.25, train_size=None),
             error_score=nan,
             estimator=KNeighborsRegressor(algorithm='auto', leaf_size=30,
                                           metric='minkowski',
                                           metric_params=None, n_jobs=None,
                                           n_neighbors=5, p=2,
                                           weights='uniform'),
             iid='deprecated', n_jobs=-1,
             param_grid=[{'metric': ['euclidean', 'minkowski'],
                          'n_neighbors': [1, 5, 10, 20],
                          'weights': ['uniform', 'distance']}],
             pre_dispatch='2*n_jobs', refit=True, return_train_score=False,
             scoring='neg_mean_squared_error', verbose=0)
In [6]:
print("Mean test score (MSE) = ",np.round(-search.cv_results_["mean_test_score"],3))
print("Std test score = ",np.round(search.cv_results_["std_test_score"],3))

Mean test score (MSE) =  [0.003 0.003 0.003 0.003 0.004 0.003 0.005 0.004 0.003 0.003 0.003 0.003
 0.004 0.003 0.005 0.004]
Std test score =  [0.002 0.002 0.001 0.001 0.002 0.002 0.002 0.002 0.002 0.002 0.001 0.001
 0.002 0.002 0.002 0.002]
In [7]:
#print(search.best_estimator_)
print("Best set of parameters : ",search.best_params_)
print("Best RMSE : ",np.round(np.sqrt(-search.best_score_),3))

Best set of parameters :  {'metric': 'euclidean', 'n_neighbors': 5, 'weights': 'distance'}
Best RMSE :  0.054

The performances of the nearest neighbor methods could be improved if the model were fit on a subset of selected variables.

In [8]:
p = X.shape[1]
rho = np.zeros(p)
for j in range(p):
    rho[j] = np.corrcoef(X[:,j],fat)[0,1]

rho10 = np.sort(rho)[p-10]
keep = np.where(rho>=rho10)[0]


search.fit(X[:,keep],fat)
Out[8]:
GridSearchCV(cv=ShuffleSplit(n_splits=100, random_state=0, test_size=0.25, train_size=None),
             error_score=nan,
             estimator=KNeighborsRegressor(algorithm='auto', leaf_size=30,
                                           metric='minkowski',
                                           metric_params=None, n_jobs=None,
                                           n_neighbors=5, p=2,
                                           weights='uniform'),
             iid='deprecated', n_jobs=-1,
             param_grid=[{'metric': ['euclidean', 'minkowski'],
                          'n_neighbors': [1, 5, 10, 20],
                          'weights': ['uniform', 'distance']}],
             pre_dispatch='2*n_jobs', refit=True, return_train_score=False,
             scoring='neg_mean_squared_error', verbose=0)
In [9]:
#print(search.best_estimator_)
print("Best set of parameters : ",search.best_params_)
print("Best RMSE : ",np.round(np.sqrt(-search.best_score_),3))

Best set of parameters :  {'metric': 'euclidean', 'n_neighbors': 1, 'weights': 'uniform'}
Best RMSE :  0.047

Regression trees

In [10]:
# Initialization of the regressor
tree = tree.DecisionTreeRegressor() 
param_grid = [{'max_features': [10,50,100]}]
# Optimization of the hyperparameters
search = GridSearchCV(tree, 
                      param_grid, 
                      scoring="neg_root_mean_squared_error",
                      cv=cv,
                      n_jobs=-1)
search.fit(X,fat)
Out[10]:
GridSearchCV(cv=ShuffleSplit(n_splits=100, random_state=0, test_size=0.25, train_size=None),
             error_score=nan,
             estimator=DecisionTreeRegressor(ccp_alpha=0.0, criterion='mse',
                                             max_depth=None, max_features=None,
                                             max_leaf_nodes=None,
                                             min_impurity_decrease=0.0,
                                             min_impurity_split=None,
                                             min_samples_leaf=1,
                                             min_samples_split=2,
                                             min_weight_fraction_leaf=0.0,
                                             presort='deprecated',
                                             random_state=None,
                                             splitter='best'),
             iid='deprecated', n_jobs=-1,
             param_grid=[{'max_features': [10, 50, 100]}],
             pre_dispatch='2*n_jobs', refit=True, return_train_score=False,
             scoring='neg_root_mean_squared_error', verbose=0)
In [11]:
#print(search.best_estimator_)
print("Best set of parameters : ",search.best_params_)
print("Best RMSE : ",-np.round(search.best_score_,3))

Best set of parameters :  {'max_features': 100}
Best RMSE :  0.052
In [12]:
res_cv = cross_val_score(tree, X, fat, scoring = "neg_root_mean_squared_error", cv=cv)
In [13]:
print(res_cv)

[-0.04278744 -0.04823902 -0.05455316 -0.06527728 -0.06335407 -0.03029454
 -0.06440289 -0.0646715  -0.0489034  -0.07176944 -0.05748207 -0.06058513
 -0.05149892 -0.0592408  -0.03985954 -0.06747626 -0.05523389 -0.04993416
 -0.04297468 -0.04043178 -0.06735035 -0.05401304 -0.05089975 -0.03780285
 -0.06329683 -0.03524367 -0.05143683 -0.04573359 -0.06107336 -0.06767015
 -0.04281357 -0.04146276 -0.05749389 -0.03955574 -0.0637104  -0.06549362
 -0.07957012 -0.03761005 -0.02715752 -0.04295817 -0.05063509 -0.05320244
 -0.04014512 -0.05991419 -0.03242303 -0.06249088 -0.05269047 -0.07355603
 -0.06375929 -0.07834459 -0.03092677 -0.07881992 -0.05360448 -0.03918389
 -0.04367491 -0.03251837 -0.04945914 -0.03609296 -0.06836985 -0.05105147
 -0.03697205 -0.04918195 -0.0633145  -0.07203387 -0.0331157  -0.04354118
 -0.05863252 -0.05372051 -0.05236538 -0.04060307 -0.04897645 -0.07504306
 -0.04547345 -0.06872299 -0.04122848 -0.06295862 -0.04299373 -0.06547654
 -0.03782754 -0.04968425 -0.0481976  -0.03566546 -0.06702455 -0.06799998
 -0.02145501 -0.05363716 -0.07729801 -0.04263078 -0.04921453 -0.06499711
 -0.07488214 -0.05065944 -0.03223709 -0.0564966  -0.07171654 -0.04857846
 -0.05463273 -0.04801898 -0.04849565 -0.04921808]
In [18]:
print("Mean RMSE for tree model = ",round(-np.mean(res_cv),3))
print("Std RMSE for tree model = ",round(np.std(res_cv),3))

Mean RMSE for tree model =  0.053
Std RMSE for tree model =  0.013
In [ ]: