# 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 _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()