import numpy as np
import math
from ..Model.estimlm import func_estimlm as estimlm
from ..basefunctions.makerow import func_makerow as makerow
from ..basefunctions.sepym import func_sepym as sepym
from ..basefunctions.sdiff import func_sdiff as sdiff
[docs]
def estimate(pmod, y, u=np.array([]), show_plot=True, show_output=True):
"""Estimate prediction-model parameters using the Levenberg–Marquardt algorithm.
Pre-processes ``y`` (and optionally ``u``) according to the transforms and
differencing orders stored in ``pmod``, then calls the LM estimator
(:func:`~TimeSeriesSRC.Model.estimlm.func_estimlm`) and returns the fitted
model together with the full training record.
Parameters
----------
pmod : pmodel
Prediction model object created by :class:`~TimeSeriesSRC.Model.model.pmodel`.
Its parameter arrays (``a``, ``b``, ``c``, ``d``, ``f``) are updated
in-place during estimation.
y : array-like, shape (1, N)
Output (response) time series. Should be zero-mean (or near zero-mean)
after applying the differencing specified in ``pmod.diff``.
u : array-like, shape (n_inputs, N), optional
Input time series matrix. Omit for purely ARMA models.
Default is an empty array.
show_plot : bool, optional
Display the live training-performance (MSE vs epoch) plot.
Default ``True``.
show_output : bool, optional
Print per-epoch training summary to stdout. Default ``True``.
Returns
-------
pmod : pmodel
The fitted model with updated parameter arrays.
trec : dict
Training record. Key ``'index'`` holds the per-epoch MSE vector;
``'epoch'`` holds the number of completed epochs.
stat : dict
Final-epoch summary statistics (MSE, gradient norm, Jacobian
condition number).
Examples
--------
>>> import numpy as np
>>> from TimeSeriesSRC.Model.model import pmodel
>>> from TimeSeriesSRC.Model.estimate import estimate
>>> y = np.array([-0.19, 0.52, -3.50, 3.01, -3.04, 1.59]).reshape(1, -1)
>>> u = np.array([-0.43, -1.67, 0.13, 0.29, -1.15, 1.19]).reshape(1, -1)
>>> pm = pmodel('bjtf', nb=[1], nc=[1], nd=[1], nf=[1], delay=[0])
>>> pm.estimParams.epochs = 5
>>> pm, trec, stat = estimate(pm, y, u, show_plot=False, show_output=False)
See Also
--------
pmodel : Prediction model object with polynomial-order specification.
func_selpmod : Automated grid search over candidate model structures.
func_pmodmse : Compute mean-squared prediction error for a fitted model.
"""
math_functions = dir(math)
pmod.set_data(y, u)
uflag = (len(u) > 0)
if uflag:
u = makerow(u)
ystru, y, m = sepym(y)
y = makerow(y)
# Preprocess the sequences
upreproc = pmod.upreproc
pr = len(upreproc)
if (uflag and pr != 0):
#if (pr != 1 and pr != u.shape[1]):
# xerror = 'rows of upreproc should either equals 1 or the number of inputs. '
# raise Exception(xerror)
for i in range(pr):
#u = eval(upreproc[pr, i], u)
if upreproc[i] in math_functions:
code = 'math.{}(x)'.format(upreproc[i])
else:
code = '{}(x)'.format(upreproc[i])
for j in range(len(u)):
uj = list(u[j])
uj = list(map(lambda x: eval(code, globals(), {'x': x}), uj))
uj = np.array(uj)
u[j] = uj
ypreproc = pmod.ypreproc
pc = len(ypreproc) # only one output is possible
#if (pc>1):
# xerror = 'ypreproc should have only one row. '
# raise Exception(xerror)
if pc != 0:
for i in range(pc):
if ypreproc[i] in math_functions:
code = 'math.{}(x)'.format(ypreproc[i])
else:
code = '{}(x)'.format(ypreproc[i])
for j in range(len(y)):
yj = list(y[j])
print(yj)
yj = list(map(lambda x: eval(code, globals(), {'x': x}), yj))
yj = np.array(yj)
y[j] = yj
# Difference the sequences before estimation so the optimizer minimises MSE
# on the stationary (differenced) series. predict() expects pre-differenced
# data; callers that need predictions on the original scale must difference
# their y (and u) before calling predict().
period = [x for x in pmod.period]
period.insert(0, 1)
diff = pmod.diff
for i in range(len(diff)):
d = diff[i]
if d != 0:
if uflag:
u = sdiff(u, d, period[i])
y = sdiff(y, d, period[i])
# check to see if y, u are zero mean
if abs(np.mean(y)) > (2 * np.std(y)):
print('The desired output may not be a zero mean sequence.');
if uflag and any(abs(np.mean(u, 1)) > (2 * np.std(u[:,1]))):
print('Input may not be zero mean sequences.');
# Call the appropriate estimation function
ystru['y'] = y
ystru['m'] = ystru['m'][0, :len(y[0])].reshape(1,-1)
if uflag:
if pmod.estimFcn == 'estimlm':
pmod, trec, stat = estimlm(pmod, ystru, u, show_plot=show_plot, show_output=show_output)
else:
if pmod.estimFcn == 'estimlm':
pmod, trec, stat = estimlm(pmod, ystru, show_plot=show_plot, show_output=show_output)
return pmod, trec, stat