Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

DOC add example for basic usage of NARX #32

Merged
merged 1 commit into from
Jan 3, 2025
Merged

DOC add example for basic usage of NARX #32

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
151 changes: 151 additions & 0 deletions examples/plot_narx.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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()
48 changes: 24 additions & 24 deletions fastcan/narx.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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
Expand All @@ -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 |
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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_)
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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,
Expand All @@ -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")],
Expand All @@ -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.
Expand All @@ -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 |
Expand Down Expand Up @@ -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
----------
Expand Down Expand Up @@ -810,8 +810,8 @@ def make_narx(

Returns
-------
narx : Narx
Narx instance.
narx : NARX
NARX instance.

Examples
--------
Expand Down Expand Up @@ -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)
Loading
Loading