Uncertainty from variable to derivative in data fitting

Question

If I have a system of differential equations, say $x' = f(x,p)$ and a set of data $(t_n,x_n)$, where $p$ is the set of constant parameters. I can then use a fitting method like least squared to obtain $p^{best fit}$ with some uncertainty $\delta p$. When using this $p^{best fit}$, it will result in some uncertainty $\delta x$ for the state variable $x$. By uncertainty, I mean $\delta p$ or $\delta x$ can be thought of as the standard deviation (perhaps of a Gaussian), but there is probably more technical ways of describing it.

My question is: how can I estimate $\delta x'$ (the uncertainty in the rate of change of $x$ with respect to time $t$)?

Thus far, I thought of these possibilities.

1. Estimate the rate data from $(t_n,x_n)$ to get $(t_m,x'_m)$ and fit the rate equation to the rate data similar to this Stack post

2. I could also use $\delta p$ to explicitly calculate the uncertainty in $x'$ with standard uncertainty propagation rules wiki. The Monte Carlo sampling method probably works as well if I assume some distributions.

3. Also, I thought using Ito calculus is another possibility to account for the randomness (in the sense of uncertainty).

Of course, these methods have their drawbacks. What is the standard approach for this kind of question?

p.s: I apologize in advance for not using proper statistical language.

Answer 1

It's not entirely clear what your procedure actually is, but I'll write this answer assuming that it's something like the following (let's assume for simplicity of notation that there are three spatial variables $x,y,z$ and three parameters $a,b,c$):

You have a system of differential equations:

$$x'=f(x,y,z,a,b,c)$$ $$y'=g(x,y,z,a,b,c)$$ $$z'=h(x,y,z,a,b,c)$$

This system of differential equations has a solution:

$$x=F(t,a,b,c)$$ $$y=G(t,a,b,c)$$ $$z=H(t,a,b,c)$$

which depends both on time and on the parameters $a,b,c$.

You also have a numerical solver that can, for a particular value of $a,b,c$, generate $F,G,H$ from $f,g,h$.

To determine the best fit, you compare $F(t)$, $G(t)$ and $H(t)$ to some input data $x_0(t)$, $y_0(t)$, and $z_0(t)$, varying $a,b,c$ to try to minimize some measure of the difference between the numerical solution and the input data. From there, you get your optimal parameters.

So what you're actually trying to do is find the uncertainty in a function of three variables and three functions of four variables: namely, you're trying to find the uncertainty in the following quantity:

$$f(x(t,a,b,c),y(t,a,b,c),z(t,a,b,c),a,b,c)$$

which has both explicit and implicit dependence on the parameters.

Generally, we cannot assume that the parameters are uncorrelated. If you're using any kind of decent fitting software, you should be able to get a covariance matrix $C$ which contains information about the correlations, defined as follows:

$$C=\begin{pmatrix}\sigma_a^2&\sigma_{ab}&\sigma_{ac}\\\sigma_{ba}&\sigma_b^2&\sigma_{cb}\\\sigma_{ca}&\sigma_{cb}&\sigma_c^2\end{pmatrix}$$

where, for example, $\sigma_a^2$ is the variance of $a$, and $\sigma_{ab}$ is the covariance between $a$ and $b$. From there, we can get the uncertainties $\sigma_x$, $\sigma_y$, and $\sigma_z$ using the solution $F,G,H$, the elements of this matrix, and the timestep $\sigma_t$:

$$\sigma_x^2=\left(\frac{\partial F}{\partial t}\right)^2\sigma_t^2+\left(\frac{\partial F}{\partial a}\right)^2\sigma_a^2+\left(\frac{\partial F}{\partial b}\right)^2\sigma_b^2+\left(\frac{\partial F}{\partial c}\right)^2\sigma_c^2+2\frac{\partial F}{\partial a}\frac{\partial F}{\partial b}\sigma_{ab}+2\frac{\partial F}{\partial a}\frac{\partial F}{\partial c}\sigma_{ac}+2\frac{\partial F}{\partial b}\frac{\partial F}{\partial c}\sigma_{bc}$$

and similarly for $y$ and $z$, substituting $G$ and $H$ respectively for $F$.

The part I'm less sure about comes in here: we also cannot assume that the uncertainties in $x,y,z$ are uncorrelated, so we need to propagate the covariances as well. I believe this is done in the following way:

$$\sigma_{xy}=\left\vert\frac{\partial F}{\partial t}\right\vert\left\vert\frac{\partial G}{\partial t}\right\vert\sigma_t^2+\left\vert\frac{\partial F}{\partial a}\right\vert\left\vert\frac{\partial G}{\partial a}\right\vert\sigma_a^2+\left\vert\frac{\partial F}{\partial b}\right\vert\left\vert\frac{\partial G}{\partial b}\right\vert\sigma_b^2+\left\vert\frac{\partial F}{\partial c}\right\vert\left\vert\frac{\partial G}{\partial c}\right\vert\sigma_c^2+\left(\frac{\partial F}{\partial a}\frac{\partial G}{\partial b}+\frac{\partial G}{\partial a}\frac{\partial F}{\partial b}\right)\sigma_{ab}+\left(\frac{\partial F}{\partial a}\frac{\partial G}{\partial c}+\frac{\partial G}{\partial a}\frac{\partial F}{\partial c}\right)\sigma_{ac}+\left(\frac{\partial F}{\partial b}\frac{\partial G}{\partial c}+\frac{\partial G}{\partial b}\frac{\partial F}{\partial c}\right)\sigma_{bc}$$

and, again, similarly for the other pairs of variables.

So now we have our variances and covariances for all of our variables in the main function, which means we can finally evaluate $\sigma_{x'}$:

$$\sigma_{x'}^2=\left(\frac{\partial f}{\partial x}\right)^2\sigma_x^2+\left(\frac{\partial f}{\partial y}\right)^2\sigma_y^2+\left(\frac{\partial f}{\partial z}\right)^2\sigma_z^2+2\frac{\partial f}{\partial x}\frac{\partial f}{\partial y}\sigma_{xy}+2\frac{\partial f}{\partial x}\frac{\partial f}{\partial z}\sigma_{xz}+2\frac{\partial f}{\partial y}\frac{\partial f}{\partial z}\sigma_{yz}+\left(\frac{\partial f}{\partial a}\right)^2\sigma_a^2+\left(\frac{\partial f}{\partial b}\right)^2\sigma_b^2+\left(\frac{\partial f}{\partial c}\right)^2\sigma_c^2+2\frac{\partial f}{\partial a}\frac{\partial f}{\partial b}\sigma_{ab}+2\frac{\partial f}{\partial a}\frac{\partial f}{\partial c}\sigma_{ac}+2\frac{\partial f}{\partial b}\frac{\partial f}{\partial c}\sigma_{bc}$$

and similarly for $\sigma_{y'}$ and $\sigma_{z'}$, substituting $g$ and $h$ for $f$. It should be fairly straightforward to see how to generalize this to higher numbers of parameters and higher dimensions of solution vector.