The GPoM package enables to explore large numbers
of models in polynomial ODEs form in order to identify
which model can reproduce the observed behaviour.
The present vignette aims to explain how to retrieve
the information obtained from the gPoMo
function.
In practice, the gPoMo
function will generate
a list of variables hereafter called outputGPoM
:
data("Ross76") # time vector tin <- Ross76[seq(1, 3000, by = 8), 1] # single time series data <- Ross76[seq(1, 3000, by = 8), 3] # global modelling # results are put in list outputGPoM outputGPoM <- gPoMo(data, tin=tin, dMax = 2, nS=c(3), show = 0, method = 'rk4', nPmin = 3, nPmax = 12, IstepMin = 400, IstepMax = 401)
To test the potential of a given model to reproduce the observed dynamics, the numerical integrability of each model is tested and the resulting behaviour may either diverge, converge to a trivial attractor (fixed points, periodic cycles) or to a chaotic attractor.
The aim of the present vignette is to show how to have
a rapid overview of the models explored or obtained
with the gPoMo
function.
In order to keep a memory of the analyzed signal, the data
used as an input are kept in $tin
(for the input time)
and $inputdata
(for the input time series).
The filtered version obtained when computing
the derivatives is put in $tfilt
(the vector time of
the filtered time series which length is lower due
to some lost at the boundaries) and in $filtdata
,
a matrix with the filtered time series (first column)
and its successive derivatives (second and next columns).
# plot data and models time series plot(outputGPoM$tin, outputGPoM$inputdata, xlab='t', ylab='y(t)', main = 'Time series', cex=0.4) lines(outputGPoM$tfiltdata, outputGPoM$filtdata[,1], type='l', col='green') legend(0,-5, c("original", "filtered"), col=c('black', 'green'), lty=c(1,NA), pch=c(NA,1), cex = 0.8)
The original differential phase portrait can thus be plotted as follows:
# plot data and models time series plot(outputGPoM$filtdata[,1], outputGPoM$filtdata[,2], xlab='y', ylab= expression(dy/dt), main = 'First projection', type = 'l', cex = 0.4) plot(outputGPoM$filtdata[,1], outputGPoM$filtdata[,3], xlab='y', ylab= expression(d^2*y/dt^2), main = 'Second projection', type = 'l', cex = 0.4)
Information about the tested models are kept in the other
variables: $okMod
, $stockoutreg
and $models
.
Variable $okMod
provides information about the obtained
numerical integration of all the tested models: values are
chosen equal to zero when the integrated trajectory
diverged during the numerical integration,
equal to 1 when the model could not be categorized
(meaning that the model is still under categorizing
test).
Since we are mainly interested in non-diverging models,
these are expected to be detected in the $okMod
parameter
such as: $okMod != 0
. Note that diverging models may be
identified as non-diverging if the numerical integration
was not performed on a sufficiently long duration.
The number of $okMod
identified as non-diverging
(based on the chosen maximum integration steps IStepMax
)
is obtained by:
# How many models could not be identified as non-diverging sum(outputGPoM$okMod)
The reference numbers of the uncategorized models can be retrieved as follows:
# what is the reference number of these models? which(outputGPoM$okMod == 1)
The numerical integration of the tested models is kept
in variable $stockoutreg
. For example, it can be plotted
for models #7 and #9:
# plot data and models time series plot(outputGPoM$tin, outputGPoM$inputdata, xlab='t', ylab='y(t)', cex = 0.4, main = 'Model 7') lines(outputGPoM$tfiltdata, outputGPoM$filtdata[,1], type='l', col='green') lines(outputGPoM$tout, outputGPoM$stockoutreg$model7[,1], type='l', col='red') legend(0,-4.5, c("original", "filtered", "model"), col=c('black', 'green', 'red'), lty=c(NA,1,1), pch=c(1, NA, NA), cex = 0.6) # plot(outputGPoM$tin, outputGPoM$inputdata, xlab='t', ylab='y(t)', cex = 0.4, main = 'Model 9') lines(outputGPoM$tfiltdata, outputGPoM$filtdata[,1], type='l', col='green') lines(outputGPoM$tout, outputGPoM$stockoutreg$model9[,1], type='l', col='red') legend(0,-4.5, c("original", "filtered", "model"), col=c('black', 'green', 'red'), lty=c(NA,1,1), pch=c(1, NA, NA), cex = 0.6)
In this example, only one iteration of the model test
was performed (since IstepMin
$* 2 <$ IstepMax
).
As a consequence, the same initial conditions
are used for both filtered data and model reconstruction:
head(outputGPoM$filtdata[1:3], 1) head(outputGPoM$stockoutreg$model7, 1)
(note that $filtdata
has one more derivative
which was necessary to identify the model
parameters).
Original and models phase portraits can also be plotted for comparison as follows:
# plot data and models phase portraits plot(outputGPoM$filtdata[,1], outputGPoM$filtdata[,2], type='l', col='green', xlab='y', ylab='dy/dt', main = 'Model 7') lines(outputGPoM$stockoutreg$model7[,1], outputGPoM$stockoutreg$model7[,2], type='l', col='red') legend(0,3.9,c("model", "filtered"), col=c('red', 'green'), lty=c(1,1), cex=0.6) # plot(outputGPoM$filtdata[,1], outputGPoM$filtdata[,2], type='l', col='green', xlab='y', ylab='dy/dt', main = 'Model 9') lines(outputGPoM$stockoutreg$model9[,1], outputGPoM$stockoutreg$model9[,2], type='l', col='red') legend(0,3.9,c("model", "filtered"), col=c('red', 'green'), lty=c(1,1), cex=0.6)
Finally, models coefficients are all kept in the sublist
$models
. These are separated in $models$mToTest#
where #
should be replaced by the model number
(the list of the models that were tested with
the gPoMo
algorithm) and $models$model#
(the list of uncategorized models).
In both cases, models coefficients are in a matrix
of dimension pMax * nVar
with pMax
the maximal
number of parameter (see vignette 1 Conventions
) and
nVar
the number of variables (also called model
dimension).
The corresponding equations can be edited directly
using the visuEq
function as:
nVar <- dim(outputGPoM$models$model7)[2] pMax <- dim(outputGPoM$models$model7)[1] dMax <- p2dMax(nVar,pMaxKnown = pMax) # See equations visuEq(outputGPoM$models$model7, approx = 2) # visuEq(outputGPoM$models$model9, approx = TRUE)
By default, a large number of digits is provided for each
parameter.
To make its reading easier, this number can be controlled
by using the optional parameter approx = TRUE
,
the minimum number of digit will be kept (the same for
all the parameters) but keeping all the terms in the
polynomial.
The number of digits to add to the minimum number of digits
required can also be chosen manually.
The number of terms of each model can be obtained as follows:
# for model #7 sum(outputGPoM$models$model7 != 0) # for model #9 sum(outputGPoM$models$model9 != 0)
Since these models have a canonical structure, it may also be interesting to know the number of terms in the last equation (which polynomial defines the model), such as:
# for model #7 colSums(outputGPoM$models$model7 != 0) # for model #9 colSums(outputGPoM$models$model9 != 0)
Finally, to have an overview of the overall outputs
of gPoMo
at a glance for all the models, the following
command can be used:
###################### # automatic ########## ###################### visuOutGP(outputGPoM)
An overview of all the main elements of the global
modelling technique and GPoM tools were presented
in the previous vignettes.
An illustration of model validation will now be
presented in the last vignette 5 Predictability
based on the model forecasting performances.
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.