I've been trying to fit some data with polynomials using a pseudo-inversion library (numpy.linalg). There are two real functions $f$ and $g$ such that $f=-\alpha g$, $\alpha > 0$, and
$$ f=f(t_1,\ldots,t_5), $$ $$ g = g(t_1,\ldots,t_5), $$
where $t_1,\ldots,t_5$ are real parameters. The parameter $\alpha$ has been estimated by an experiment, so I have two sets of experimental values for $f$ and $g$, $\{f_1,\ldots,f_{12}\}$, $\{g_1,\ldots,g_{12}\}$.
What I want to do is, given a good model, to reproduce this experimental result ($\alpha$) numerically.
Starting from rather a good prior estimation of the parameters, I generate five set of numbers for each parameter:
- for the first one, $\{(t_1)_1, (t_1)_2,\ldots, (t_1)_{12}\}$;
- for the second one, $\{(t_2)_1, (t_2)_2,\ldots, (t_2)_{12}\}$;
- etc...
so that $(t_i)_{j+k} - (t_i)_{j} \equiv (t_a)_{b+k} - (t_a)_{b}$ $\forall i,j,a,b,k$, that is varying them by the same differential every time.
Next, I write the problem to first order like this
$$ \left[ \begin{matrix} f_1\\ f_2\\ \vdots\\ f_{12} \end{matrix} \right] = \left[ \begin{matrix} 1& (t_1)_1 & (t_2)_1 & \dots & (t_5)_1 \\ 1& (t_1)_2 & (t_2)_2 & \dots & (t_5)_2 \\ \vdots\\ 1& (t_1)_{12} & (t_2)_{12} & \dots & (t_5)_{12} \\ \end{matrix} \right] \left[ \begin{matrix} c_0\\ c_1\\ \vdots\\ c_5 \end{matrix} \right] $$
or
$$ \vec{f} = T \cdot\vec{c}. $$
I then pseudo-invert $T$ and apply the result to $f$ and get a good estimation for $\vec{c} = [c_0,c_1,\ldots]$, that I will call $\vec{c}_f$.
I repeat for $g$ and get $\vec{c}_g$.
Of course $\vec c_g\neq \vec c_f$, however
$$ \frac{\vec c_f}{c_{f1}} \simeq \frac{\vec c_g}{c_{g1}}, $$
which represents
$$ \frac{f}{\partial f/\partial t_1} \simeq \frac{c_{f0}}{\partial f/\partial t_1} + t_1 + \frac{\partial t_1}{\partial t_2} t_2 + \dots $$
or
$$ \frac{f}{c_{f1}} \simeq \frac{c_{f0}}{c_{f1}} + t_1 + \frac{c_{f2}}{c_{f1}} t_2 + \dots $$
with this I know that, to first order,
$$ dt_2 \simeq dt_1 \times \frac{c_{f1}}{c_{f2}} $$
so I can finally run my model where, for each unit of variation $dt_1$ of the first parameter, the $n^{th}$ one varies this much:
$$ dt_n \simeq dt_1 \times \frac{c_{f1}}{c_{fn}}, $$
and the result is nice, but I can do better.
Question number one: how to I set up this calculation (pseudoinversion) to higher order? I presume that having higher order coefficients is sufficient, but please confirm, I'm insecure.
Question number two: It's not easy to build this matrix (higher order), so are there python libraries that you are aware of that build such matrix automatically given a dataset? If not, I should generate every possible combination of mixed derivatives I think, right?
Question number three (the most important): Let's say I do construct the correct matrix (not easy) and find the best fit coefficients,what do I do with them?
To second order, for instance, I have something like (please remember that all of this is finite)
$$ \frac{f}{\partial f/\partial t_1} \simeq \frac{c_{f0}}{\partial f/\partial t_1} + t_1 + \frac{\partial t_1}{\partial t_2} t_2 + \frac{1}{2} \frac{\partial t_1}{\partial f} \frac{\partial^2 f}{\partial t_2\partial t_2} t_2^2 + \dots $$
where
$$ \frac{\partial t_1}{\partial f} \frac{\partial^2 f}{\partial t_2\partial t_2} = \frac{\partial}{\partial t_2} \frac{\partial t_1}{\partial t_2} - \frac{\partial f}{\partial t_2} \frac{\partial }{\partial f}\frac{\partial t_1}{\partial t_2} $$
and similar equations can be written for mixed derivatives terms.
How do I use these coefficients to determine how $t_n$ varies with respect to $t_1$, that is writing an equation similar to this:
$$ dt_n \simeq dt_1 \times \frac{c_{f1}}{c_{fn}}, $$
but for higher orders?
Thank you for reading and thank you very much in advance for your help
Edit: My shy and inconclusive attempt at second order corrections
If I define
$$ c_{\alpha_1\dots\alpha_p} = \frac{\Delta(\Delta\dots(\Delta f)\dots)}{\Delta t_{\alpha_1}\dots \Delta t_{\alpha_p}} $$
I get
$$ \Delta t_\beta = \frac{c_\alpha}{c_\beta} \Delta t_\alpha $$
for first order and
$$ \Delta t_\beta = \frac{c_{\gamma\delta}}{c_{\alpha\beta}} \frac{\Delta t_\gamma \Delta t_\delta}{\Delta t_\alpha} $$
for second order terms.
This implies
$$ \Delta t_\gamma \Delta t_\delta = \frac{c_\beta c_{\alpha\beta}}{c_\alpha c_{\gamma\delta}} $$
which one may consider a second order term, but is it?
Please correct me if I'm wrong with what follows: to solve the first order problem I had to fix a normalization term and I chose $\alpha = 1$ by dividing the Taylor expansion by
$$ \partial f/\partial t_1. $$
Now that I'm dealing with higher order, I'm apparently required to fix two normalization terms instead, so by choosing
$$ \alpha = \beta = 1, $$
I get, for $\delta = \gamma$,
$$ (\Delta t_\gamma)^2 = \frac{c_{11}}{c_{\gamma\gamma}} $$
that confuses me, since it's in no way related to, for instance, $\Delta t_1$ or any other finite differential.
Is it ok? Should I expect that the correct way to modulate any of this parameters could be, to second order, independent of what I do with any other?
Is it even right to call this a second order correction? How would I even modulate $t_\gamma$ given these two (first and second order) contributions? Just by adding them together with a $1/2!$ factor in front of the latter?
All of this sounds wrong to me but I have no idea how to tackle this.
How?
Thank you again