DOC add example for basic of NARX (#32)

Author: MatthewSZhang

examples/plot_narx.py

@@ -0,0 +1,151 @@
+"""
+===============================================
+Nonlinear AutoRegressive eXogenous (NARX) model
+===============================================
+
+.. currentmodule:: fastcan
+
+In this example, we illustrate how to build a polynomial
+NARX model for time series prediction.
+"""
+
+# Authors: Sikai Zhang
+# SPDX-License-Identifier: MIT
+
+# %%
+# Prepare data
+# ------------
+#
+# First, a simulated time series dataset is generated from the following nonlinear
+# system.
+#
+# .. math::
+#     y(k) = 0.5y(k-1) + 0.3u_0(k)^2 + 2u_0(k-1)u_0(k-3) + 1.5u_0(k-2)u_1(k-3) + 1
+#
+# where :math:`k` is the time index,
+# :math:`u_0` and :math:`u_1` are input signals,
+# and :math:`y` is the output signal.
+
+
+import numpy as np
+
+rng = np.random.default_rng(12345)
+n_samples = 1000
+max_delay = 3
+e = rng.normal(0, 0.1, n_samples)
+u0 = rng.uniform(0, 1, n_samples+max_delay)
+u1 = rng.normal(0, 0.1, n_samples+max_delay)
+y = np.zeros(n_samples+max_delay)
+for i in range(max_delay, n_samples+max_delay):
+    y[i] = 0.5*y[i-1]+0.3*u0[i]**2+2*u0[i-1]*u0[i-3]+1.5*u0[i-2]*u1[i-3]+1
+y = y[max_delay:]+e
+X = np.c_[u0[max_delay:], u1[max_delay:]]
+
+# %%
+# Build term libriary
+# -------------------
+#
+# To build a polynomial NARX model, it is normally have two steps:
+#
+# 1. Search the structure of the model, i.e., the terms in the model, e.g.,
+# :math:`u_0(k-1)u_0(k-3)`, :math:`u_0(k-2)u_1(k-3)`, etc.
+# 2. Learn the coefficients before the terms.
+#
+# To search the structure of the model, the candidate term libriary should be
+# constructed.
+#
+# 1. Time-shifted variables: the raw input-output data, i.e., :math:`u0(k)`,
+# :math:`u1(k)`, and :math:`y(k)`, are converted into :math:`u0(k-1)`, :math:`u1(k-2)`,
+# etc.
+# 2. Nonlinear terms: the time-shifted variables are onverted to nonlinear terms
+# via polynomial basis functions, e.g., :math:`u_0(k-1)^2`, :math:`u_0(k-1)u_0(k-3)`,
+# etc.
+#
+#   .. rubric:: References
+#
+#   * `"Nonlinear system identification: NARMAX methods in the time, frequency,
+#     and spatio-temporal domains" <https://doi.org/10.1002/9781118535561>`_
+#     Billings, S. A. John Wiley & Sons, (2013).
+#
+#
+# Make time-shifted variables
+# ^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+from fastcan.narx import make_time_shift_features, make_time_shift_ids
+
+time_shift_ids = make_time_shift_ids(
+    n_features=3, # Number of inputs (2) and output (1) signals
+    max_delay=3, # Maximum time delays
+    include_zero_delay = [True, True, False], # Whether to include zero delay
+    # for each signal. The output signal should not have zero delay.
+)
+
+time_shift_vars = make_time_shift_features(np.c_[X, y], time_shift_ids)
+
+# %%
+# Make nonlinear terms
+# ^^^^^^^^^^^^^^^^^^^^
+
+from fastcan.narx import make_poly_features, make_poly_ids
+
+poly_ids = make_poly_ids(
+    n_features=time_shift_vars.shape[1], # Number of time-shifted variables
+    degree=2, # Maximum polynomial degree
+)
+
+poly_terms = make_poly_features(time_shift_vars, poly_ids)
+
+# %%
+# Term selection
+# --------------
+# After the term library is constructed, the terms can be selected by :class:`FastCan`,
+# whose :math:`X` is the nonlinear terms and :math:`y` is the output signal.
+
+from fastcan import FastCan
+
+selector = FastCan(
+    n_features_to_select=4, # 4 terms should be selected
+).fit(poly_terms, y)
+
+support = selector.get_support()
+selected_poly_ids = poly_ids[support]
+
+
+# %%
+# Build NARX model
+# ----------------
+# As the polynomical NARX is a linear function of the nonlinear tems,
+# the coefficient of each term can be easily estimated by oridnary least squares.
+
+from fastcan.narx import NARX, print_narx
+
+narx_model = NARX(
+    time_shift_ids=time_shift_ids,
+    poly_ids = selected_poly_ids,
+)
+
+narx_model.fit(X, y)
+
+# In the printed NARX model, it is found that :class:`FastCan` selects the correct
+# terms and the coefficients are close to the true values.
+print_narx(narx_model)
+
+# %%
+# Plot NARX prediction performance
+# --------------------------------
+
+import matplotlib.pyplot as plt
+from sklearn.metrics import r2_score
+
+y_pred = narx_model.predict(
+    X[:100],
+    y_init=y[:narx_model.max_delay_] # Set the initial values of the prediction to
+    # the true values
+)
+
+plt.plot(y[:100], label="True")
+plt.plot(y_pred, label="Predicted")
+plt.xlabel("Time index k")
+plt.legend()
+plt.title(f"NARX prediction R-squared: {r2_score(y[:100], y_pred):.5f}")
+plt.show()

fastcan/narx.py

@@ -274,9 +274,9 @@ def _mask_missing_value(*arr):
     return tuple([x[mask_nomissing] for x in arr])
 
 
-class Narx(RegressorMixin, BaseEstimator):
+class NARX(RegressorMixin, BaseEstimator):
     """The Nonlinear Autoregressive eXogenous (NARX) model class.
-    For example, a (polynomial) Narx model is like
+    For example, a (polynomial) NARX model is like
     y(t) = y(t-1)*u(t-1) + u(t-1)^2 + u(t-2) + 1.5
     where y(t) is the system output at time t,
     u(t) is the system input at time t,
@@ -333,7 +333,7 @@ class Narx(RegressorMixin, BaseEstimator):
     Examples
     --------
     >>> import numpy as np
-    >>> from fastcan.narx import Narx, print_narx
+    >>> from fastcan.narx import NARX, print_narx
     >>> rng = np.random.default_rng(12345)
     >>> n_samples = 1000
     >>> max_delay = 3
@@ -351,7 +351,7 @@ class Narx(RegressorMixin, BaseEstimator):
     >>> poly_ids = [[0, 2], # 1*u(k-1)
     ...             [0, 4], # 1*y(k-1)
     ...             [1, 3]] # u(k-1)*u(k-3)
-    >>> narx = Narx(time_shift_ids=time_shift_ids,
+    >>> narx = NARX(time_shift_ids=time_shift_ids,
     ...             poly_ids=poly_ids).fit(X, y, coef_init="one_step_ahead")
     >>> print_narx(narx)
     |        Term        |   Coef   |
@@ -397,16 +397,16 @@ def fit(self, X, y, sample_weight=None, coef_init=None, **params):
 
         sample_weight : array-like of shape (n_samples,), default=None
             Individual weights for each sample, which are used for a One-Step-Ahead
-            Narx.
+            NARX.
 
         coef_init : array-like of shape (n_terms,), default=None
             The initial values of coefficients and intercept for optimization.
-            When `coef_init` is None, the model will be a One-Step-Ahead Narx.
+            When `coef_init` is None, the model will be a One-Step-Ahead NARX.
             When `coef_init` is `one_step_ahead`, the model will be a Multi-Step-Ahead
-            Narx whose coefficients and intercept are initialized by the a
-            One-Step-Ahead Narx.
+            NARX whose coefficients and intercept are initialized by the a
+            One-Step-Ahead NARX.
             When `coef_init` is an array, the model will be a Multi-Step-Ahead
-            Narx whose coefficients and intercept are initialized by the array.
+            NARX whose coefficients and intercept are initialized by the array.
 
             .. note::
                 When coef_init is None, missing values (i.e., np.nan) are allowed.
@@ -477,7 +477,7 @@ def fit(self, X, y, sample_weight=None, coef_init=None, **params):
         n_terms = self.poly_ids_.shape[0] + 1
 
         if isinstance(coef_init, (type(None), str)):
-            # fit a one-step-ahead Narx model
+            # fit a one-step-ahead NARX model
             xy_hstack = np.c_[X, y]
             osa_narx = LinearRegression()
             time_shift_vars = make_time_shift_features(xy_hstack, self.time_shift_ids_)
@@ -508,7 +508,7 @@ def fit(self, X, y, sample_weight=None, coef_init=None, **params):
                 )
 
         lsq = least_squares(
-            Narx._residual,
+            NARX._residual,
             x0=coef_init,
             args=(
                 self._expression,
@@ -578,7 +578,7 @@ def _residual(
         coef = coef_intercept[:-1]
         intercept = coef_intercept[-1]
 
-        y_hat = Narx._predict(expression, X, y[:max_delay], coef, intercept, max_delay)
+        y_hat = NARX._predict(expression, X, y[:max_delay], coef, intercept, max_delay)
 
         y_masked, y_hat_masked = _mask_missing_value(y, y_hat)
 
@@ -622,7 +622,7 @@ def predict(self, X, y_init=None):
                     f"but got {y_init.shape}."
                 )
 
-        return Narx._predict(
+        return NARX._predict(
             self._expression,
             X,
             y_init,
@@ -639,7 +639,7 @@ def __sklearn_tags__(self):
 
 @validate_params(
     {
-        "narx": [Narx],
+        "narx": [NARX],
         "term_space": [Interval(Integral, 1, None, closed="left")],
         "coef_space": [Interval(Integral, 1, None, closed="left")],
         "float_precision": [Interval(Integral, 0, None, closed="left")],
@@ -652,12 +652,12 @@ def print_narx(
     coef_space=10,
     float_precision=3,
 ):
-    """Print a Narx model as a Table which contains Term and Coef.
+    """Print a NARX model as a Table which contains Term and Coef.
 
     Parameters
     ----------
-    narx : Narx model
-        The Narx model to be printed.
+    narx : NARX model
+        The NARX model to be printed.
 
     term_space: int, default=20
         The space for the column of Term.
@@ -671,14 +671,14 @@ def print_narx(
     Returns
     -------
     table : str
-        The table of terms and coefficients of the Narx model.
+        The table of terms and coefficients of the NARX model.
 
     Examples
     --------
     >>> from sklearn.datasets import load_diabetes
-    >>> from fastcan.narx import print_narx, Narx
+    >>> from fastcan.narx import print_narx, NARX
     >>> X, y = load_diabetes(return_X_y=True)
-    >>> print_narx(Narx().fit(X, y), term_space=10, coef_space=5, float_precision=0)
+    >>> print_narx(NARX().fit(X, y), term_space=10, coef_space=5, float_precision=0)
     |   Term   |Coef |
     ==================
     |Intercept | 152 |
@@ -765,7 +765,7 @@ def make_narx(
     refine_max_iter=None,
     **params,
 ):
-    """Find `time_shift_ids` and `poly_ids` for a Narx model.
+    """Find `time_shift_ids` and `poly_ids` for a NARX model.
 
     Parameters
     ----------
@@ -810,8 +810,8 @@ def make_narx(
 
     Returns
     -------
-    narx : Narx
-        Narx instance.
+    narx : NARX
+        NARX instance.
 
     Examples
     --------
@@ -912,4 +912,4 @@ def make_narx(
         - 1
     )
 
-    return Narx(time_shift_ids=time_shift_ids, poly_ids=poly_ids)
+    return NARX(time_shift_ids=time_shift_ids, poly_ids=poly_ids)