Source code for sysidentpy.residues.residues_correlation

# Authors:
#           Wilson Rocha Lacerda Junior <wilsonrljr@outlook.com>
#           Luan Pascoal da Costa Andrade <luan_pascoal13@hotmail.com>
#           Samuel Carlos Pessoa Oliveira <samuelcpoliveira@gmail.com>
#           Samir Angelo Milani Martins <martins@ufsj.edu.br>
# License: BSD 3 clause


import numpy as np
import matplotlib.pyplot as plt


[docs]class ResiduesAnalysis: """Residues analysis for Polynomial NARX model."""
[docs] def residuals(self, X, y, yhat): """Performs the residual analysis of output to validate model. Parameters ---------- y : array-like of shape = n_samples The target data used in the identification process. yhat : array-like of shape = n_samples The prediction values of the identification process. X : ndarray of floats The input data. Returns ------- output_autocorr : ndarray of floats: 1st column - Residuals normalized autocorrelation. 2nd/3rd columns - Superior and inferior limits of a 95% confidence interval. output_crosscorr : ndarray of floats: 1st column - Correlation between residuals and input. 2nd/3rd columns - Superior and inferior limits of a 95% confidence interval. References ---------- [1] `Wikipedia entry on the Autocorrelation <https://en.wikipedia.org/wiki/Autocorrelation>`_ Examples -------- >>> y = [3, -0.5, 2, 7] >>> autocorr(y) [62.25 11.5 2.5 21. ] """ e, e_acf, unnormalized_e_acf = self._residuals_acf(y, yhat) xe_ccf = self._input_ccf(X[:, 0], e, len(unnormalized_e_acf)) ex = e * X[:, 0] ee = e**2 x2line = (X[:, 0]**2) - (X[:, 0]**2).mean() eeline = (e**2) - (e**2).mean() e_ex = self._input_ccf(e, ex, len(unnormalized_e_acf)) x2line_e = self._input_ccf(x2line, e, len(unnormalized_e_acf)) x2line_ee = self._input_ccf(x2line, ee, len(unnormalized_e_acf)) ye = y * e.reshape(-1, 1) yeline = (y * e.reshape(-1, 1)) - (y * e.reshape(-1, 1)).mean() yeline_x2line = self._input_ccf(yeline, x2line, len(unnormalized_e_acf)) lam = np.sqrt(sum((ee - ee.mean())**2) / sum((ye - ye.mean())**2)) yeline_eeline = self._input_ccf(yeline, eeline, len(unnormalized_e_acf)) return (e_acf, xe_ccf, [yeline_x2line, yeline_eeline, e_ex, x2line_e, x2line_ee], lam)
[docs] def _input_ccf(self, X, e, len_confidence): xe_ccf = np.zeros((int(np.floor(len(X)/2)), 3), dtype=float) xe_ccf[:, 0] = self._normalized_correlation(X, e) xe_ccf[:, 1] = np.ones( len(xe_ccf[:, 0])) * (1.96/np.sqrt(len_confidence)) xe_ccf[:, 2] = xe_ccf[:, 1] * (-1) return xe_ccf
[docs] def _residuals_acf(self, y, yhat): e = (y - yhat).flatten() unnormalized_e_acf = np.correlate(e, e, mode='full') half_of_simmetry_autocorr = int(np.floor(unnormalized_e_acf.size/2)) e_acf = np.zeros( (len(unnormalized_e_acf)-half_of_simmetry_autocorr, 3), dtype=float) e_acf[:, 0] = unnormalized_e_acf[half_of_simmetry_autocorr:] / \ unnormalized_e_acf[half_of_simmetry_autocorr] e_acf[:, 1] = np.ones( len(e_acf[:, 0])) * (1.96/np.sqrt(len(unnormalized_e_acf))) e_acf[:, 2] = e_acf[:, 1] * (-1) return e, e_acf, unnormalized_e_acf
[docs] def _normalized_correlation(self, signal1, signal2): """Compute the normalized correlation between two signals. Parameters ---------- signal1 : array-like of shape = n_samples. signal2 : array-like of shape = n_samples. Returns ------- ruy : ndarray of floats: The normalized cross correlation between the two signals. """ y = (signal1 - np.mean(signal1)).flatten() u = (signal2 - np.mean(signal2)).flatten() t = int(np.floor(len(signal1)/2)) ruy = np.array(np.zeros(t)) ruy[0] = np.sum(y*u) / (np.sqrt(np.sum(y**2)) * np.sqrt(np.sum(u**2))) for i in range(1, t): y = (signal1 - np.mean(signal1[:i])).flatten() u = (signal2 - np.mean(signal2[i:])).flatten() ruy[i] = np.sum( y[:-i]*u[i:]) / (np.sqrt(np.sum(y[:-i]**2)) * np.sqrt(np.sum(u[i:]**2))) return ruy
[docs] def plot_result(self, y, yhat, e_acf, xe_ccf, figsize=(10, 8), n=100): """Plot the free run simulation and residues analysis. Parameters ---------- y : array-like of shape = n_samples The target data used in the identification process. yhat : array-like of shape = n_samples The prediction values of the identification process. e_acf : ndarray of floats: 1st column - Residuals normalized autocorrelation. 2nd/3rd columns - Superior and inferior limits of a 95% confidence interval. xe_ccf : ndarray of floats: 1st column - Correlation between residuals and input. 2nd/3rd columns - Superior and inferior limits of a 95% confidence interval. """ plt.style.use('seaborn-whitegrid') plt.rcParams['axes.facecolor'] = 'white' fig, axes = plt.subplots(nrows=2, ncols=2, figsize=figsize, facecolor='white') fig.subplots_adjust(hspace=0.7) for ax, feature in zip(axes.flatten()[2:], [e_acf, xe_ccf]): ax.plot(feature) ax.set_xlabel('Lag', fontsize=12) ax.set_ylabel('Cross Correlation: ee, ex', fontsize=12) # ax = plt.gca() ax.set_ylim([-1, 1]) ax.grid(color='grey', linestyle='-.', alpha=0.2) ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) ax.xaxis.set_ticks_position('bottom') ax.yaxis.set_ticks_position('left') fig.tight_layout() ax = plt.subplot(211) ax.plot(y[self.max_lag:n], c='#101820FF', alpha=0.2, marker='o', label='Data', linewidth=5) ax.plot(yhat[self.max_lag:n], c='#006B38FF', marker='*', label='Model', linewidth=1.) ax.set_title('Free run simulation', fontsize=18) ax.legend() ax.tick_params(labelsize=14) ax.set_xlabel('Samples', fontsize=14) ax.set_ylabel('y, yhat', fontsize=14) ax.grid(color='grey', linestyle='-.', alpha=0.2) ax.spines['right'].set_visible(False) ax.spines['top'].set_visible(False) ax.xaxis.set_ticks_position('bottom') ax.yaxis.set_ticks_position('left') fig.tight_layout() plt.show()