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.

  • $\begingroup$ Just to be clear: you are not going to vary the $c_j$ constants, right? If so you do not need to account for the variation of $m_i$ with them. $\endgroup$
    – Alexis
    Jan 25 '19 at 22:21
  • $\begingroup$ The $c_j$ are are essentially fixed values, coming from a database. But they have some error, which I must not neglect. $\endgroup$
    – jitter
    Jan 26 '19 at 11:04

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.

  • $\begingroup$ Thank you. This approach sounds promising, since it avoids any hassle with the derivatives naturally. My model function $m$ is also not insanely expensive, so MCMC should be applicable. I will give it a try. Do you have perhaps any reading recommendation, particularly with practical hints for an MCMC user? $\endgroup$
    – jitter
    Jan 26 '19 at 11:17
  • $\begingroup$ An overview can be found here. I personally would recommend the book from Andrew Gelman, which is also listed here. $\endgroup$
    – Bort
    Jan 27 '19 at 11:35

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.