gloMoId {GPoM}

Global Model Identification

Description

Algorithm for global modelling in polynomial and canonical formulation of Ordinary Differential Equations. Univariate Global modeling aims to obtain multidimensional models from single time series (Gouesbet & Letellier 1994, Mangiarotti et al. 2012). An example of such application can be found in Mangiarotti et al. (2014) For a multivariate application, see GPoMo (Mangiarotti 2015, Mangiarotti et al. 2016).

Example:
For a model dimension nVar=3, the global model will read:
dX1/dt = X2
dX2/dt = X3
dX3/dt = P(X1,X2,X3).

Usage

gloMoId(
  series,
  tin = NULL,
  dt = NULL,
  nVar = NULL,
  dMax = 1,
  dMin = 0,
  weight = NULL,
  show = 1,
  filterReg = NULL,
  winL = 9
)

Arguments

series

The original data set: either a single vector corresponding to the original variable; Or a matrix containing the original variable in the first column and its successive derivatives in the next columns. In the latter case, for the construction of n-dimensional model, series should have nVar+1 columns since one more derivative will be necessary to identify the model parameters. Variable nVar will be set equal to n. In the former case, that is when only a single vector is provided, the derivatives will be automatically recomputed. Therefore, the dimension nVar expected for the model has to be provided.

tin

Input date vector which length should correspond to the input time series.

dt

Sampling time of the input time series.

nVar

Number of variables considered in the polynomial formulation.

dMax

Maximum degree of the polynomial formulation.

dMin

The minimum negative degree of the polynomial formulation (0 by default).

weight

A vector providing the binary weighting function of the input data series (0 or 1). By default, all the values are set to 1.

show

Provide (2) or not (0-1) visual output during the running process.

filterReg

A vector that specifies the template for the equation structure (for one single equation). The convention defined by poLabs is used. Value is 1 if the regressor is available, 0 if it is not.

winL

Total number of points used for computing the derivatives of the input time series. This parameter will be used as an input in function drvSucc to compute the derivatives.

Value

A list of five elements :

$init The original time series and the successive derivatives used for the modeling.

$filterReg The structure of the output model. Value is 1 if the regressor is available, 0 if it is not. The terms order is given by function poLabs.

$K Values of the identified coefficients corresponding to the regressors defined in filterReg.

$resTot The variance of the residual signal of the model.

$resSsMod The variance of the residual signal of the closer submodels.

$finalWeight Weighting series after boundary values were removed.

Author(s)

Sylvain Mangiarotti, Laurent Drapeau, Mireille Huc

References

[1] Gouesbet G., Letellier C., Global vector-field reconstruction by using a multivariate polynomial L2 approximation on nets, Physical Review E, 49 (6), 4955-4972, 1994.
[2] Mangiarotti S., Coudret R., Drapeau L., & Jarlan L., Polynomial search and global modeling : Two algorithms for modeling chaos, Physical Review E, 86, 046205, 2012.
[3] Mangiarotti S., Drapeau L. & Letellier C., Two chaotic models for cereal crops observed from satellite in northern Morocco. Chaos, 24(2), 023130, 2014.
[4] Mangiarotti S., Low dimensional chaotic models for the plague epidemic in Bombay (1896-1911), Chaos, Solitons & Fractals, 81(A), 184-196, 2015.
[5] Mangiarotti S., Peyre M. & Huc M., A chaotic model for the epidemic of Ebola Virus Disease in West Africa (2013-2016). Chaos, 26, 113112, 2016.

See Also

gPoMo, autoGPoMoSearch, autoGPoMoTest, poLabs

Examples

#############
# Example 1 #
#############
# load data
data("Ross76")
tin <- Ross76[,1]
data <- Ross76[,2:3]

# Polynomial identification
reg <- gloMoId(data[0:500,2], dt=1/100, nVar=2, dMax=2, show=0)

#############
# Example 2 #
#############
# load data
data(NDVI)

# Definition of the Model structure
terms <- c(1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1)
poLabs(3,3)[terms==1]
reg <- gloMoId(NDVI [,1:1], dt=1/125, nVar=3, dMax=3,
               show=0, filterReg=terms==1)



#############
# Example 3 #
#############
# load data
data("Ross76")
# time vector
tin <- Ross76[1:500,1]
# single time series
series <- Ross76[1:500,3]
# some noise is added
series[1:100] <- series[1:100] + 0.01 * runif(1:100, min = -1, max = 1)
series[301:320] <- series[301:320] + 0.05 * runif(1:20, min = -1, max = 1)
# weighting function
W <- tin * 0 + 1
W[1:100] <- 0  # the first hundred values will not be considered
W[301:320] <- 0  # twenty other values will not be considered either
reg <- gloMoId(series, dt=1/100, weight = W, nVar=3, dMax=2, show=1)
visuEq(reg$K, 3, 2, approx = 4)
# first weight which value not equal to zero:
i1 = which(reg$finalWeight == 1)[1]
v0 <-  reg$init[i1,1:3]

reconstr <- numicano(nVar=3, dMax=2, Istep=5000, onestep=1/250, PolyTerms=reg$K,
                     v0=v0, method="ode45")
plot(reconstr$reconstr[,2], reconstr$reconstr[,3], type='l', lwd = 3,
                            main='phase portrait', xlab='time t', ylab = 'x(t)', col='orange')
# original data:
lines(reg$init[,1], reg$init[,2], type='l',
      main='phase portrait', xlab='x', ylab = 'dx/dt', col='black')
# initial condition
lines(v0[1], v0[2], type = 'p', col = 'red')