Importância do Extended Least Squares¶
Exemplo criado por Wilson Rocha Lacerda Junior
Procurando mais detalhes sobre modelos NARMAX? Para informações completas sobre modelos, métodos e uma ampla variedade de exemplos e benchmarks implementados no SysIdentPy, confira nosso livro: Nonlinear System Identification and Forecasting: Theory and Practice With SysIdentPy
Este livro oferece orientação aprofundada para apoiar seu trabalho com o SysIdentPy.
Aqui importamos o modelo NARMAX, a métrica para avaliação do modelo e os métodos para gerar dados de amostra para testes. Também importamos pandas para uso específico.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sysidentpy.model_structure_selection import FROLS
from sysidentpy.basis_function import Polynomial
from sysidentpy.parameter_estimation import LeastSquares
from sysidentpy.metrics import root_relative_squared_error
from sysidentpy.utils.generate_data import get_siso_data
from sysidentpy.utils.display_results import results
Gerando dados de amostra com 1 entrada e 1 saída¶
Os dados são gerados simulando o seguinte modelo: \(y_k = 0.2y_{k-1} + 0.1y_{k-1}x_{k-1} + 0.9x_{k-2} + e_{k}\)
Se colored_noise for definido como True:
\(e_{k} = 0.8\nu_{k-1} + \nu_{k}\)
onde \(x\) é uma variável aleatória uniformemente distribuída e \(\nu\) é uma variável com distribuição gaussiana com \(\mu=0\) e \(\sigma\) definido pelo usuário.
No próximo exemplo, geraremos dados com 3000 amostras com ruído colorido e selecionando 90% dos dados para treinar o modelo.
x_train, x_valid, y_train, y_valid = get_siso_data(
n=1000, colored_noise=True, sigma=0.2, train_percentage=90
)
Construindo o modelo¶
Primeiro, treinaremos um modelo sem o Algoritmo Extended Least Squares para fins de comparação.
basis_function = Polynomial(degree=2)
estimator = LeastSquares(unbiased=False)
model = FROLS(
order_selection=False,
n_terms=3,
ylag=2,
xlag=2,
info_criteria="aic",
estimator=estimator,
basis_function=basis_function,
err_tol=None,
)
model.fit(X=x_train, y=y_train)
yhat = model.predict(X=x_valid, y=y_valid)
rrse = root_relative_squared_error(y_valid, yhat)
print(rrse)
0.5499799245432233
Claramente, há algo errado com o modelo obtido. Veja o notebook basic_steps para comparar os resultados obtidos usando os mesmos dados mas sem ruído colorido. Mas vamos ver o que está errado.
r = pd.DataFrame(
results(
model.final_model,
model.theta,
model.err,
model.n_terms,
err_precision=8,
dtype="sci",
),
columns=["Regressors", "Parameters", "ERR"],
)
print(r)
Regressors Parameters ERR
0 x1(k-2) 8.9976E-01 7.41682256E-01
1 y(k-1) 2.8734E-01 8.33321202E-02
2 x1(k-1)y(k-1) 1.2348E-01 5.10334067E-03
Estimação de parâmetros com bias¶
Como podemos observar acima, a estrutura do modelo é exatamente a mesma que gerou os dados. Você pode ver que o ERR ordenou os termos da forma correta. E esta é uma nota importante sobre o algoritmo Error Reduction Ratio usado aqui: ele é muito robusto ao ruído colorido!!
Esta é uma grande característica! No entanto, embora a estrutura esteja correta, os parâmetros do modelo não estão ok! Aqui temos uma estimação com bias! O parâmetro real para \(y_{k-1}\) é \(0.2\), não \(0.3\).
Neste caso, estamos na verdade modelando usando um modelo NARX, não NARMAX. A parte MA existe para permitir uma estimação não-enviesada dos parâmetros. Para alcançar uma estimação não-enviesada dos parâmetros, temos o algoritmo Extended Least Squares. Lembre-se, se os dados têm apenas ruído branco, NARX é suficiente.
Antes de aplicar o Algoritmo Extended Least Squares, executaremos vários modelos NARX para verificar quão diferentes os parâmetros estimados são dos reais.
parameters = np.zeros([3, 50])
for i in range(50):
x_train, x_valid, y_train, y_valid = get_siso_data(
n=3000, colored_noise=True, train_percentage=90
)
model.fit(X=x_train, y=y_train)
parameters[:, i] = model.theta.flatten()
# Definir o tema para seaborn (opcional)
sns.set_theme()
plt.figure(figsize=(14, 4))
# Plotar KDE para cada parâmetro
sns.kdeplot(parameters.T[:, 0], label="Parameter 1")
sns.kdeplot(parameters.T[:, 1], label="Parameter 2")
sns.kdeplot(parameters.T[:, 2], label="Parameter 3")
# Plotar linhas verticais onde os valores reais devem estar
plt.axvline(x=0.1, color="k", linestyle="--", label="Real Value 0.1")
plt.axvline(x=0.2, color="k", linestyle="--", label="Real Value 0.2")
plt.axvline(x=0.9, color="k", linestyle="--", label="Real Value 0.9")
plt.xlabel("Parameter Value")
plt.ylabel("Density")
plt.title("Kernel Density Estimate of Parameters")
plt.legend()
plt.show()
Usando o algoritmo Extended Least Squares¶
Como mostrado na figura acima, temos um problema para estimar o parâmetro para \(y_{k-1}\). Agora usaremos o Algoritmo Extended Least Squares.
No SysIdentPy, basta definir extended_least_squares como True e o algoritmo será aplicado.
basis_function = Polynomial(degree=2)
estimator = LeastSquares(unbiased=True)
parameters = np.zeros([3, 50])
for i in range(50):
x_train, x_valid, y_train, y_valid = get_siso_data(
n=3000, colored_noise=True, train_percentage=90
)
model = FROLS(
order_selection=False,
n_terms=3,
ylag=2,
xlag=2,
elag=2,
info_criteria="aic",
estimator=estimator,
basis_function=basis_function,
)
model.fit(X=x_train, y=y_train)
parameters[:, i] = model.theta.flatten()
plt.figure(figsize=(14, 4))
# Plotar KDE para cada parâmetro
sns.kdeplot(parameters.T[:, 0], label="Parameter 1")
sns.kdeplot(parameters.T[:, 1], label="Parameter 2")
sns.kdeplot(parameters.T[:, 2], label="Parameter 3")
# Plotar linhas verticais onde os valores reais devem estar
plt.axvline(x=0.1, color="k", linestyle="--", label="Real Value 0.1")
plt.axvline(x=0.2, color="k", linestyle="--", label="Real Value 0.2")
plt.axvline(x=0.9, color="k", linestyle="--", label="Real Value 0.9")
plt.xlabel("Parameter Value")
plt.ylabel("Density")
plt.title("Kernel Density Estimate of Parameters")
plt.legend()
plt.show()
Ótimo! Agora temos uma estimação não-enviesada dos parâmetros!
Nota¶
Nota: O Extended Least Squares é um algoritmo iterativo. No SysIdentPy, o padrão é 30 iterações (uiter=30) porque é conhecido da literatura que o algoritmo converge rapidamente (cerca de 10 ou 20 iterações).