gPoMo: Generalized Polynomial Modeling

Description Usage Arguments Value Author(s) References See Also Examples

Description

Algorithm for a Generalized Polynomial formulation of multivariate Global Modeling. Global modeling aims to obtain multidimensional models from single time series [1-2]. In the generalized (polynomial) formulation provided in this function, it can also be applied to multivariate time series [3-4].

Example:
Note that nS provides the number of dimensions used from each variable

case I
For nS=c(2,3) means that 2 dimensions are reconstructed from variable 1: the original variable X1 and its first derivative X2), and 3 dimensions are reconstructed from variable 2: the original variable X3 and its first and second derivatives X4 and X5. The generalized model will thus be such as:
dX1/dt = X2
dX2/dt = P1(X1,X2,X3,X4,X5)
dX3/dt = X4
dX4/dt = X5
dX5/dt = P2(X1,X2,X3,X4,X5).

case II
For nS=c(1,1,1,1) means that only the original variables X1, X2, X3 and X4 will be used. The generalized model will thus be such as:
dX1/dt = P1(X1,X2,X3,X4)
dX2/dt = P2(X1,X2,X3,X4)
dX3/dt = P3(X1,X2,X3,X4)
dX4/dt = P4(X1,X2,X3,X4).

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
gPoMo(
  data,
  tin = NULL,
  dtFixe = NULL,
  dMax = 2,
  nS = c(3),
  winL = 9,
  weight = NULL,
  show = 1,
  verbose = 1,
  underSamp = NULL,
  EqS = NULL,
  AndManda = NULL,
  OrMandaPerEq = NULL,
  IstepMin = 2,
  IstepMax = 2000,
  nPmin = 1,
  nPmax = 14,
  tooFarThr = 4,
  FxPtThr = 1e-08,
  LimCyclThr = 1e-06,
  nPminPerEq = 1,
  nPmaxPerEq = NULL,
  method = "rk4"
)

Arguments

data

Input Time series: Each column is one time series that corresponds to one variable.

tin

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

dtFixe

Time step used for the analysis. It should correspond to the sampling time of the input data. Note that for very large and very small time steps, alternative units may be used in order to stabilize the numerical computation.

dMax

Maximum degree of the polynomial formulation.

nS

A vector providing the number of dimensions used for each input variables (see Examples 1 and 2). The dimension of the resulting model will be nVar = sum(nS).

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.

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.

verbose

Gives information (if set to 1) about the algorithm progress and keeps silent if set to 0.

underSamp

Number of points used for undersampling the data. For undersamp = 1 the complete time series is used. For undersamp = 2, only one data out of two is kept, etc.

EqS

Model template including all allowed regressors. Each column corresponds to one equation. Each line corresponds to one polynomial term as defined by function poLabs.

AndManda

AND-mandatory terms in the equations (all the provided terms should be in the equations).

OrMandaPerEq

OR-mandatory terms per equations (at least one of the provided terms should be in each equation).

IstepMin

The minimum number of integration step to start of the analysis (by default IstepMin = 10).

IstepMax

The maximum number of integration steps for stopping the analysis (by default IstepMax = 10000).

nPmin

Corresponds to the minimum number of parameters (and thus of polynomial term) allowed.

nPmax

Corresponds to the maximum number of parameters (and thus of polynomial) allowed.

tooFarThr

Divergence threshold, maximum value of the model trajectory compared to the data standard deviation. By default a trjactory is too far if the distance to the center is larger than four times the variance of the input data.

FxPtThr

Threshold used to detect fixed points.

LimCyclThr

Threshold used to detect the limit cycle.

nPminPerEq

Corresponds to the minimum number of parameters (and thus of polynomial term) allowed per equation.

nPmaxPerEq

Corresponds to the maximum number of parameters (and thus of polynomial) allowed per equation.

method

The integration technique used for the numerical integration. By default, the fourth-order Runge-Kutta method (method = 'rk4') is used. Other methods such as 'ode45' or 'lsoda' may also be chosen. See package deSolve for details.

Value

A list containing:

$tin The time vector of the input time series

$inputdata The input time series

$tfiltdata The time vector of the filtered time series (boudary removed)

$filtdata A matrix of the filtered time series with its derivatives

$okMod A vector classifying the models: diverging models (0), periodic models of period-1 (-1), unclassified models (1).

$coeff A matrix with the coefficients of one selected model

$models A list of all the models to be tested $mToTest1, $mToTest2, etc. and all selected models $model1, $model2, etc.

$tout The time vector of the output time series (vector length corresponding to the longest numerical integration duration)

$stockoutreg A list of matrices with the integrated trajectories (variable X1 in column 1, X2 in 2, etc.) of all the models $model1, $model2, etc.

Author(s)

Sylvain Mangiarotti, Flavie Le Jean, Mireille Huc

References

[1] Gouesbet G. & Letellier C., 1994. Global vector-field reconstruction by using a multivariate polynomial L2 approximation on nets, Physical Review E, 49 (6), 4955-4972.
[2] Mangiarotti S., Coudret R., Drapeau L. & Jarlan L., Polynomial search and Global modelling: two algorithms for modeling chaos. Physical Review E, 86(4), 046205.
[3] Mangiarotti S., Le Jean F., Huc M. & Letellier C., Global Modeling of aggregated and associated chaotic dynamics. Chaos, Solitons and Fractals, 83, 82-96.
[4] S. Mangiarotti, M. Peyre & M. Huc, 2016. A chaotic model for the epidemic of Ebola virus disease in West Africa (2013-2016). Chaos, 26, 113112.

See Also

gloMoId, autoGPoMoSearch, autoGPoMoTest

autoGPoMoSearch, autoGPoMoTest, visuOutGP, poLabs, predictab, drvSucc

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
#Example 1
data("Ross76")
tin <- Ross76[,1]
data <- Ross76[,3]
dev.new()
out1 <- gPoMo(data, tin = tin, dMax = 2, nS=c(3), show = 1,
              IstepMax = 1000, nPmin = 9, nPmax = 11)
visuEq(out1$models$model1, approx = 4)


## Not run: 
#Example 2
data("Ross76")
tin <- Ross76[,1]
data <- Ross76[,3]
# if some data are not valid (vector 'weight' with zeros)
W <- tin * 0 + 1
W[1:100] <- 0
W[700:1500] <- 0
W[2000:2800] <- 0
W[3000:3500] <- 0
dev.new()
out2 <- gPoMo(data, tin = tin, weight = W,
                 dMax = 2, nS=c(3), show = 1,
                 IstepMax = 6000, nPmin = 9, nPmax = 11)
visuEq(out2$models$model3, approx = 4)

## End(Not run)


## Not run: 
#Example 3
data("Ross76")
tin <- Ross76[,1]
data <- Ross76[,2:4]
dev.new()
out3 <- gPoMo(data, tin=tin, dMax = 2, nS=c(1,1,1), show = 1,
              IstepMin = 10, IstepMax = 3000, nPmin = 7, nPmax = 8)
# the simplest model able to reproduce the observed dynamics is model #5
visuEq(out3$models$model5, approx = 3, substit = 1) # the original Rossler system is thus retrieved

## End(Not run)

## Not run: 
#Example 4
data("Ross76")
tin <- Ross76[,1]
data <- Ross76[,2:3]
# model template:
EqS <- matrix(1, ncol = 3, nrow = 10)
EqS[,1] <- c(0,0,0,1,0,0,0,0,0,0)
EqS[,2] <- c(1,1,0,1,0,1,1,1,1,1)
EqS[,3] <- c(0,1,0,0,0,0,1,1,0,0)
visuEq(EqS, substit = c('X','Y','Z'))
dev.new()
out4 <- gPoMo(data, tin=tin, dMax = 2, nS=c(2,1), show = 1,
      EqS = EqS, IstepMin = 10, IstepMax = 2000,
      nPmin = 9, nPmax = 11)
visuEq(out4$models$model2, approx = 2, substit = c("Y","Y2","Z"))

## End(Not run)

## Not run: 
#Example 5
# load data
data("TSallMod_nVar3_dMax2")
#multiple (six) time series
tin <- TSallMod_nVar3_dMax2$SprK$reconstr[1:400,1]
TSRo76 <- TSallMod_nVar3_dMax2$R76$reconstr[,2:4]
TSSprK <- TSallMod_nVar3_dMax2$SprK$reconstr[,2:4]
data <- cbind(TSRo76,TSSprK)[1:400,]
dev.new()
# generalized Polynomial modelling
out5 <- gPoMo(data, tin = tin, dMax = 2, nS = c(1,1,1,1,1,1),
              show = 0, method = 'rk4',
              IstepMin = 2, IstepMax = 3,
              nPmin = 13, nPmax = 13)

# the original Rossler (variables x, y and z) and Sprott (variables u, v and w)
# systems are retrieved:
visuEq(out5$models$model347, approx = 4,
       substit = c('x', 'y', 'z', 'u', 'v', 'w'))
# to check the robustness of the model, the integration duration
# should be chosen longer (at least IstepMax = 4000)

## End(Not run)

GPoM documentation built on Feb. 18, 2020, 5:08 p.m.