Robustly estimating exponential growth?

Question

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.

Answer 1

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,

$$c_1e^{-c_2}+z$$

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$).

$$c_1=\dfrac{\displaystyle\sum_{k=1}^ne^{-c_2x_k}f_k}{\displaystyle\sum_{k=1}^ne^{-2c_2x_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.

Answer 2

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.

Answer 3

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.

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