Say we have a model $$f(x) = c_1e^{-c_2x}$$ We know we have some noise, but we don't know where it is, if inside exponential or added outside and we don't know the distribution.

How can we robustly estimate $c_1,c_2$ from a sequence of measured data points $$\{(x_1,f(x_1)+e_1),\cdots,(x_n,f(x_n)+e_n)\}$$ Here we denote the additive error $e_n$ even though we don't know if the noise is additive or not.

  • $\begingroup$ What do you mean by robustly in this context ? $\endgroup$ – Yves Daoust Jul 1 at 21:12
  • $\begingroup$ The approach should not fail too badly in the presence of ( possibly non gaussian ) noise. $\endgroup$ – mathreadler Jul 2 at 0:02
  • $\begingroup$ Usually, robust means in the presence of outliers, which is a different thing. You are probably asking for an unbiased solution. $\endgroup$ – Yves Daoust Jul 2 at 6:35
  • $\begingroup$ Not necessarily. Outliers can be members of some distribution too. $\endgroup$ – mathreadler Jul 2 at 7:15
  • 1
    $\begingroup$ Of course, outliers have a distribution. Not necessarily what ? $\endgroup$ – Yves Daoust Jul 2 at 7:17

Without more information, let us assume that the noise (let $z$) is homoscedastic and zero-mean.

If the model is

$$c_1e^{-c_2x+z},$$ taking the logarithm linearizes to

$$\log c_1-c_2 x+z,$$ which is easily solved by ordinary, linear least-squares.

Now if the noise is purely additive,


you have to minimize

$$E:=\sum_{k=1}^n \left(c_1e^{-c_2x_k}-f_k\right)^2,$$

leading to the nasty system of equations

$$\begin{cases}\displaystyle\sum_{k=1}^n e^{-c_2x_k}\left(c_1e^{-c_2x_k}-f_k\right)=0, \\\displaystyle\sum_{k=1}^n x_ke^{-c_2x_k}\left(c_1e^{-c_2x_k}-f_k\right)=0.\end{cases}$$

A possible approach is to let $c_2$ vary, and to draw $c_1$ from the first equation (as a weighted average of the $f_k$).


This way, you can evaluate $E$ as a function of $c_2$, and use a 1D minimizer, or use a 1D equation solver on the second equation. The $c_2$ value obtained from the other model is probably a good initial value.


You could state the different models

$$(1)\qquad y = a\exp(bx)+\varepsilon$$ $$(2)\qquad y = a\exp(bx+\varepsilon)$$ $$(3)\qquad y = a\exp(bx\varepsilon)$$ $$(4)\qquad y = a\exp(bx)\varepsilon.$$

Now, we solve for the error $\varepsilon$ to obtain

$$(1)\qquad \varepsilon = y - a\exp(bx)$$ $$(2)\qquad \varepsilon = \ln \dfrac y a - bx$$ $$(3)\qquad \varepsilon = \dfrac 1 {bx} \ln \dfrac y a$$ $$(4)\qquad \varepsilon = \dfrac y a \exp(-bx).$$

Now, use the loss function

$$J(a,b) = \sum_{n=1}^{N}\varepsilon^2_n + \lambda_1\left[a^2+b^2\right]+\lambda_2\left[|a|+|b|\right]$$

which is a standard squared loss function with elastic net regularization. Split your dataset into a training set and test data set. Then apply a hyperparameter optimization on the training set to determine the optimal values for $\lambda_1$ and $\lambda_2$ using k-fold cross-validation. Finally, evaluate the models by using the test data set and choose the model that has the lowest loss.

  • 1
    $\begingroup$ Your models 3 and 4 cannot work with a zero-mean noise. $\endgroup$ – Yves Daoust Jul 2 at 7:16

My own overly complicated approach for this problem would be like this:

We solve an operator equation $$ \bf Pd = 0 \Leftrightarrow \min_P \|Pd\|_2^2$$

$\bf d$ is our known data vector, $\bf P$ is our unknown operator.

We need to impose extra constraints regularizations because the problem is in general very underdetermined.

So we require that $\bf P$ be time-invariant, depend only on one time step and only as a multiplicative way.

Like this:

$$\bf P = \begin{bmatrix}\bf p&\bf -I&0&0&\cdots\\0&\bf p&\bf -I&0&\cdots\\\vdots&0&\ddots&\ddots&\ddots\\0&\cdots&0&\bf p&\bf -I\end{bmatrix}$$

For block-matrices $\bf p$ and $\bf I$ being identity. So this data point multiplied by "something" ($\bf p$) equals next data point.

And as you probably guessed at this point our block matrices for this particular problem are chosen to be the famous $2\times 2$ matrix representation for complex number multiplication:

$$z = a+bi \text{ represented by } \begin{bmatrix} a&-b\\b&a \end{bmatrix}$$

Plot to show that it works : Blue is data point and red cross is predicted by model.

enter image description here

Here is zoomed in upper left corner of regularized solution $\bf P$:

enter image description here


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.