Correct weighting in least squares fitting

Question

I am trying to fit some data points $d_i$ to a non-linear model function $m_i$, which depends on a number of fit parameters $f_k$ (I want to determine these) and also on some known, constant values $c_j$. For this I am using the Levenberg–Marquardt (LM) algorithm, which iteratively minimizes

\begin{equation} \chi^2=\sum_i\left[\frac{d_i-m_i\left(f_1,f_2,\ldots,c_1,c_2,\ldots\right)}{\sigma_i\left(d_i\right)}\right]^2\,,\tag{1} \end{equation}

where $\sigma_i^2\left(d_i\right)$ is the variance of $d_i$ (I'll call $\sigma_i$ the error of $d_i$). I have an estimate for the errors, so let's say they are known.

Now, I want to deal with the "constant" values $c_j$, which actually come from different measurements and thus have some (known) error themselves. My first idea is to adjust the term in the denominator, with the simplest type of error propagation:

\begin{equation} \chi^2=\sum_i\left[\frac{d_i-m_i\left(f_1,f_2,\ldots,c_1,c_2,\ldots\right)}{\sqrt{\sigma_i^2\left(d_i\right)+\sum_j\left(\frac{\partial{}m_i}{\partial{}c_j}\sigma\left(c_j\right)\right)^2}}\right]^2\,.\tag{2} \end{equation}

My question is how to calculate the partial derivatives $\frac{\partial{}m_i}{\partial{}c_j}$ correctly, since they depend on the current value of the fit parameters $f_k$. I see two solutions, but I don't find either very satisfying:

1) Apply the LM algorithm multiples times. For example: Calculate the derivatives with a starting guess $f_k^0$ and execute the LM fit, giving the solution $f_k^1$. Then, use $f_k^1$ to calculate the derivatives and execute a second LM fit.

2) Adjust the LM algorithm code to actually recalculate the denominator after each LM-iteration with the current best estimates of $f_k$.

How would you recommend I approach this problem? Any hints in general or reading suggestions are also appreciated.

Answer 1

So the correct way is to use the probability distribution of your $c$ in the modeling process. Imho, the most natural way of describing and solving this type of problem is the

Bayesian approach

$$P(f|X,m)=\frac{P(X|f,m)P(f|m)}{P(X|m)}.$$ Noise free $c$

You have the data points $X_i=(x_i,d_i)$. The model $m$ predicts the values of $d$: $$d(x)=m(x,f)+\epsilon$$ with parameters $f$ and e.g. $\epsilon =d-m=\mathcal{N}(0,\sigma)$. You assume $x$ is noise free but $d$ is not. The likelihood for observing one data point $X_i$ is $$P(X_i|f,m)=\frac{1}{\sigma\sqrt{2\pi}}e^{-\frac{\left(d_i-m(x_i,f)\right)^2}{2\sigma^2}}$$

and the log likelihood for all data points is $$\ln(P(X|f,m))=\sum_i \ln\left(P(d_i|x_i,f,m)\right),$$ which leads to your minimization expression $$\ln(P(X|f,m)) \propto \sum_i\left[\frac{d_i-m_i\left(f_1,f_2,\ldots,c_1,c_2,\ldots\right)}{\sigma}\right]^2.$$

For the parameters $f$ you are interested in the posterior pdf $$P(f|X,m)= \frac{P(X|f,m)P(f|m)}{P(X|m)}.$$ This can be solved by assuming some prior for $f$ and doing markov chain monte carlo (MCMC) sampling.

The model evidence $P(X|m)$ is the likelihood averaged over the model parameters and tells how well the data are predicted given the model. $$P(X|m)=\int P(X,f|m)df=\int P(X|f,m)P(f|m)df$$

Noisy $c$

Now you do the same exercise again with the additional parameter $c$ and construct the posterior pdf for $f$ given $c$, $X$, $m$, also the parameter distribution $P(c)=\mathcal{N}(c_m,\sigma_c)$ is known and independent of $f$ and $m$.

The posterior pdf $P(f|X,m)$ is then the marginalization of $P(f|X,c,m)$ $$P(f|X,m)=\int P(f|X,c,m)P(c)dc.$$

In Bayesian modeling there are several frameworks which provide MCMC sampling also for hierarchical models. e.g. PyMC is capable of doing the MCMC calculation together with an additional distribution for $c$. There is however a drawback. If the model $m$ is very expensive, sampling can quickly get out of hands.