I'm working on a nonlinear multivariate autoregressive model of order 1 (markovian). It is a discrete-time dynamical system which models exchange of mass between compartments in a compartmental model where each compartment corresponds a given population type characterized, for every time instant, by a fraction of mass $x_u(t)$. The network dynamics are expressed by quantities of mass which have chances of interact and being converted to another type, thus leading to coupled interactions which determine flow of mass between types.

The mass for type $u$ at time $t + 1$ is given by its mass at $t$ plus the incoming flows, coming from a subset of neighboring types, minus outgoing flow to other types.

$$ x_u(t+1) = x_u(t) + \sum_{v : u \in N_{\mathcal C}(v)} f_{vu}(t) -\sum_{v \in \mathcal N_{\mathcal C}(u)} f_{uv}(t) $$

Where each flow has the form:

$$ f_{uv}(t) = x_u (t) \cdot \sum_{z \in S_u} \frac{x_z(t)}{(\alpha_{u} - 1 )\cdot x_u(t) + B } \cdot \frac{p_{uz}}{m} $$

My objective is doing estimation of parameters $\alpha_u$ (1 for each compartment) and $p_{uv}$ (one for each couple of compartments) from a time series of observations concerning the masses $x_u$. I obtained good results using least squares nonlinear parameters estimation for low dimensional systems, while it quickly becomes impractical for systems with moderate to high number of compartments. Could someone suggest a technique/direction to learn the parameters of this kind of models?

EDIT: the parameters $p_{uv}$ are constrained to be probabilities (thus ranging in [0,1]) and $\alpha_u \geq 1$


Let's denote $\hat{x}_u(t)$ the mass predicted by the model given the parameter estimates $\hat{\alpha}_u$ and $\hat{p}_{uv}$. Then, you can define the following function:

$$ S(\alpha,p_{uv};x_u) = \sum_{u,t} \left( x_u(t) - \hat{x}_u(t) \right)^2 \,. $$

This function returns a scalar that is the sum over all $u$ and $t$ of the squared differences between the observed values and the predicted values. Thus, it is a measure of the prediction errors given a set of parameter estimates. (The square is taken to avoid that negative errors are offset positive errors.)

Then, the goal is to minimize this error function. This could be done by means of a general purpose optimization algorithm, for example: the Nelder-Mead algorithm or the Broyden–Fletcher–Goldfarb–Shanno algorithm. Starting from an arbitrary set of parameters, these algorithms follow an iterative procedure seeking for the values of the parameters that minimises the objective function. The values of the parameters found by this approach will be optimal in the sense that minimise the prediction errors.

Note: If there is more than one local minima, there is no guarantee that these algorithms will converge to a global minimum. Running the algorithms from different starting values is advisable.

  • $\begingroup$ I already tried with PEM nonlinear optimization suite from MATLAB (using the Trust-Region-Reflective algorithm for parameter space exploration) which minimizes the prediction error (as per this question). Again, the problem is the high dimensionality which requires this kind of objective functions (with n(n+1) parameters) to be infeasible. Do you think Nelder-Mead or Broyden–Fletcher–Goldfarb–Shanno are more computationally convenient? $\endgroup$ – maru Jan 27 '15 at 8:27
  • $\begingroup$ I don't know how these algorithms would work in your case. The point of my answer was to define a general framework and strategy to learn the parameters of the model. The R package BB contains some functions for the optimization of large-scale nonlinear systems. The algorithms and references used there may give you further ideas. Adding the tag optimization may result in more views and potential answers. $\endgroup$ – javlacalle Jan 27 '15 at 17:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.