3

Assuming I have the following model:

$$ y(t) = \alpha {e}^{- \beta t} + \gamma + n(t) $$ Where $ n(t) $ is additive white Gaussian noise (AWGN) and $ \alpha, \beta, \gamma $ are the unknown parameters to be estimated.

In this case linearization using the logarithm function won't help.

What could I do? What would be the ML Estimator / LS Method? Is there a special treatment to non positive data points?

How about the following model:

$$ y(t) = \alpha {e}^{- \beta t} + n(t) $$

In this case, it would be helpful to use the logarithm function, yet, I could I handle the non positive data?

Thanks.

cardinal
  • 24,973
  • 8
  • 94
  • 128
Royi
  • 966
  • 8
  • 20
  • 1
    It helps to define initialisms like AWGN (*additive white Gaussian noise*, I presume) upon first use, as they tend to be particular to certain fields or subfields. – cardinal Dec 27 '11 at 13:00
  • Have you looked in, e.g., the classical [Seber & Wild](http://www.amazon.com/Nonlinear-Regression-Wiley-Probability-Statistics/dp/0471471356)? – cardinal Dec 27 '11 at 13:01
  • Some general discussion of models like these appears at http://stats.stackexchange.com/questions/6720/estimation-of-exponential-model. – whuber Dec 27 '11 at 15:23
  • @cardinal, I will look into the book. Could you point me specifically to this problem? Thanks. – Royi Dec 28 '11 at 09:34

1 Answers1

4

I don’t really get what you mean with AGWN, is this simply that $n(t_i)$ are independent with $n(t_i) \sim N(0,\sigma^2)$?

The least squares estimator (which, as usual with a normal model, coincide with maximum likelihood estimator, see eg this answer) is easy to find with numerical methods, here is a piece of R code:

a <- 2; b <- 0.2; c <- -1
t <- seq(-5,10, length=100)
y <- a*exp(-b*t)+c+rnorm(length(t), sd=2)
f <- function(par, t, y) { 
  a <- par[1]; b <- par[2]; c <- par[3]; 
  return(sum((a*exp(-b*t)+c-y)**2 ) ); }
nlm( function(par) f(par, t, y), c(0,0,0))

The result of the last call is 

$minimum
[1] 430.8242

$estimate
[1]  2.0875336  0.1961210 -0.8672079

$gradient
[1]  1.301591e-05 -2.745537e-05  2.603429e-05

$code
[1] 1

$iterations
[1] 19

Our initial values (2, 0.2, -1) are estimated by (2.088, 0.196, -0.867). You can get a plot of the data points, the "true model" (red line) and the estimated model (dotted red line) as follows:

plot(t,y)
lines( t, a*exp(-b*t)+c, col="red")
nlm( function(par) f(par, t, y), c(0,0,0))$estimate -> r
lines( t, r[1]*exp(-r[2]*t)+r[3], col="red", lty=2)

Data points, with in red line the true model, in dotted red line the fitted model

Also, be sure not to miss this related question already quoted by whuber in the comments.

A few hints for numerical optimization

In the above code the numerical optimization is done using nlm, a function of R. Here are a few hints for a more elementary solution.

If $b\ne 0$ is fixed, it is easy to find $a$ and $c$ minimizing $f(a,b,c) = \sum_{i=1}^n (a \exp(-b\cdot t_i) +c -y_i)^2$: this a the classical least squares for the linear model. Namely, letting $x_i = \exp(-b\cdot t_i)$, we find $$ \begin{array}{rcccl} a &=& a(b) &=& { n \sum_i y_i x_i - \left( \sum_i y_i \right)\left(\sum_i x_i \right)\over n \sum x_i^2 - \left(\sum_i x_i \right)^2 }, \\ c &=& c(b) &=& { \sum_i y_i \sum_i x_i^2 - \left( \sum_i y_i x_i \right)\left(\sum_i x_i \right)\over n \sum x_i^2 - \left(\sum_i x_i \right)^2 }. \\ \end{array}$$

Then set $$ g(b) = f(a(b), b, c(b)).$$ The optimization problem is now reduced to find the minimum of $g(b)$. This can be done either using Newton method or a ternary search.

Community
  • 1
Elvis
  • 11,870
  • 36
  • 56
  • +1 It would help to explain why this really is a least squares problem. – whuber Dec 27 '11 at 15:23
  • @whuber I don’t get it, do you mean "why the MLE and the LSE coincide"? – Elvis Dec 27 '11 at 15:30
  • Yes: the steps implicit in your answer are that (a) you propose using ML and (b) the MLE is the least squares solution. With your experience, these are natural and perhaps seem obvious, but other readers might not make the connections. – whuber Dec 27 '11 at 16:28
  • 1
    You’re right @whuber and this is not so obvious, just "usual". I has the feeling that the OP was familiar with this, I added a link to a very similar question, I hope this is enough — if you think that it would be better to rewrite it in this special case let me know. – Elvis Dec 27 '11 at 17:43
  • @ElvisJaggerAbdul-Jabbar, Thank you. Could you explain the algorithm as I'm not familiar with 'R'. How could you handle "Negative" samples and the constant. – Royi Dec 28 '11 at 09:26
  • @Drazick I just minimize the function $(a,b,c) \mapsto \sum_i (a \exp(-b\cdot t_i) + c - y_i)^2$. This is done here using `nlm` which is a numerical procedure. You will find such procedures in every programming languages, if not as part of the language in function libraries. So there is nothing special with negative $y_i$. – Elvis Dec 28 '11 at 15:57
  • @Elvis, I knew I could always use LS by minimizing the cost function using grid search (And optimization to make it quicker). My question was, is there an analytic method? Some kind of quick linearization? I don't want to use exhaustive methods. Thanks. – Royi Dec 28 '11 at 18:40
  • @Drazick You are right, an exhaustive grid search would be a terrible solution. Here the `nlm` function uses a Newton like method. I had a quick look at your problem, I don’t think there is an analytic solution, however a Newton like method converges very fast (you see on the example above that the `nlm` function used "only" 19 iterations to converge; this could be surely improved by a fine tuning of the method, but it is not bad already). Please tell us what is your favorite programming language. But if you wish so I can write you some pseudo code and some R code for this particular problem. – Elvis Dec 28 '11 at 19:40
  • Alternatively, you can remark that when you fix $b$, letting $x_i = e^{-bt_i}$ leads to a linear model $y_i = a x_i + c + \epsilon_i$, and you find easily closed forms for $a(b)$ and $c(b)$, depending on $b$ (and $y_i$, $t_i$). So you just have to find $b$ for which the value $f(a(b),b,c(b))$ is minimal, this is a one dimensional optimization problem, it should be easier to deal with. – Elvis Dec 28 '11 at 21:14
  • @Elvis, I mainly use MATLAB. I don't want to use any internal function. Could you assist with that? Thank You. – Royi Dec 29 '11 at 08:12
  • @Drazick Matlab has many optimization tools eg `fminsearch`: http://www.mathworks.fr/help/techdoc/ref/fminsearch.html which looks perfect. If for some reason you don’t want to use it (but why on earth would you refuse to use your tools? I’m not sure I got what you mean by "I don't want to use any internal function", if you use Matlab you use its internal functions), use my previous comment to write a single parameter function $g(b) = f(a(b),b,c(b))$ and reprogram either Newton method for $g$, or a [ternary search](http://en.wikipedia.org/wiki/Ternary_search). I’ll add a few lines in my answer. – Elvis Dec 29 '11 at 09:28
  • @Elvis, I will mark your answer later since I am using smartphone and this function doesn't work. Anyhow, I can't use MATLAB function as I wanna be able to easily implement it in hardware. Hence I need a step by step method or analytical method. Thank You! – Royi Dec 30 '11 at 16:53
  • @Drazick I understand better now. Well I think Newton approach is ok, you should converge in few iterations; however you should take the time to compute literary $g'(b)$ and $g''(b)$ instead of using a numerical approximation; that’s tedious but the quality of the result will be improved. – Elvis Dec 31 '11 at 14:22
  • @Elvis, I will give it a deeper look next week. If I have any more question I'll come back. Thank you for your great assistance. – Royi Jan 05 '12 at 12:42