4
$\begingroup$

I am trying to estimate parameters, 4 of them, by fitting a system of 3 ordinary differential equations. I am using a model published that was using 3 parameters and gave a good fit to the data, and I am trying to add a parameter to this (without changing the system of ODE, there is just a rate that is used twice, and I decided to make it two parameters instead).

More precisely, the system in question is: $$\left\{\begin{array} yy_1' = -\alpha(k_1+k_2)y_1 + \alpha k_2y_2\\ y_2' = \beta k_1y_1 - (\beta k_2 + k_4)y_3\\ y_3' = \gamma k_4y_2 - \delta y_3 \end{array} \right.$$ In the original model, $k_1=k_2$, but I wanted to try with different values.

Since with this change, I have now more parameters to fit (the $k_i, \ i=1,...,4$) than I have dependent variables ($y_i,\ i = 1,...,3$), I was wondering if it could have been the reason why I could not get a proper fit using nonlinear least squares, in the same way the authors had done.

I was trying to minimize the sum of the squared errors, namely: $$SSE(\theta) = \frac{1}{3N}\sum_{i=1}^3\sum_{j=0}^{N} (y_i(t_j,\theta)-\bar y_i(t_j))^2$$ with $\theta=(k_1,k_2,k_3,k_4)$ the parameters and $\bar y_i(t_j)$ the observations.

Turned out it was just a problem of the nonlinear least square problem being to hard for the algorithms and that using the log of the model and observations made it better conditioned. Namely, the problem of minimizing: $$SSE_{log}(\theta) = \frac{1}{3N}\sum_{i=1}^3\sum_{j=0}^{N} \Bigl(\ln\bigl(y_i(t_j,\theta)\bigr)-\ln\bigl(\bar y_i(t_j)\bigr)\Bigr)^2$$ was a much simpler task and I obtained a satisfying fit to my data.

Improve this question
$\endgroup$
9
  • 2
    $\begingroup$ I think that a system of 3 linear ODE (of first order) can have more than 4 free parameters. $\endgroup$
    – nicoguaro ♦
    Aug 28, 2021 at 17:37
  • 1
    $\begingroup$ Could you write the system of equations that you're fitting? $\endgroup$
    – nicoguaro ♦
    Aug 28, 2021 at 18:11
  • 1
    $\begingroup$ The number of parameters should match the number of constraints which is very different from the number of ODEs in this problem. For example, if those three ODEs define a trajectory in 3D space - x(t), y(t), z(t) - and this trajectory should pass close to a given point r0=(x0,y0,z0) that is one constraint only. If the trajectory should also pass close to another point r1=(x1,y1,z1) that's your second constraint. Or it can be a cumulative constraint, like passing through a cloud of points $r_k$. But either way, the number of those ODEs is unrelated to the number of constraints. $\endgroup$
    – Maxim Umansky
    Aug 28, 2021 at 22:32
  • 4
    $\begingroup$ The question is not so much how many ODEs you have, but how many observations you have. If you have observed the solution of the ODE at 1000 points in time, then you have 3000 data points and that should be enough to identify 4 parameters, assuming the problem is well posed. $\endgroup$
    – Wolfgang Bangerth
    Aug 28, 2021 at 22:46
  • 3
    $\begingroup$ Nonlinear least squares is often a very nonlinear optimization problem and it might be difficult to find a solution. I would start with the three parameter problem and fit those parameters. Then you make that fourth parameter independent from its "twin" and you start the optimization routine again with four parameters, but using the three-parameter solution as a starting guess. $\endgroup$
    – Wolfgang Bangerth
    Aug 28, 2021 at 22:48

1 Answer 1

Reset to default
1
$\begingroup$

This seems to be a particular case of the general nonlinear inverse problem, a huge branch of Geophysics (or Statistics, you could say). There is an infinitude of methods to attack it; for example: under certain conditions, the system you mentioned may be suitable for Monte-Carlo methods, among others. It's important to understand the system - for example, "granulometry" is critical. Some physical systems (as in many electromagnetic experiments) can show singularities and several local minima in several dimensions. If the grid isn't properly chosen and based on the characteristics dimensions of the phenomena (rather of the data), you can easily miss the general optimum - and never find it. My advice would be:

  1. Get a theoretical/empirical idea about the phenomena being modeled. What is it? In particular, what are the characteristics lengths (time, space, ...)? Then you get an idea about the proper smoothness of the problem => range/sampling density for the calculations. In working with real data (fitting), remember that these choices cannot be completely independent of the data (Nyquist theorem).
  2. Gain an understanding of the system working with projections: solving sub-system for some parameters, while others remain constant under certain values (as you make some tests, you'll find interesting values for some parameters that you can set as constant to test the system behavior.
  3. If you are confident enough about your system's smoothness on the chosen grid (range and sampling), you can choose among many nonlinear inversion techniques. There is a plethora of them. Some can perform only if the system is factorable,...

For most cases, especially for dense grids (or large datasets), Monte-Carlo methods are worth try.

https://www.researchgate.net/publication/231129060_Monte_Carlo_analysis_of_inverse_problems

http://www.ipgp.fr/~tarantola/Files/Professional/Papers_PDF/MonteCarlo_latex.pdf

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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