2

I have the following model:

$y(t) = C + A_{0}\exp^{(-\sigma t)}\cos(\omega_dt+\phi)$

which I've coded in R as:

function(t,C,Ao,s,wd,ph) C + Ao * exp(-s*t) * cos(wd*t + ph)

I want to use this equation to form a predictive model.

However, I can't figure out how to either successfully run or plot this equation.

  • I tried nls but got various errors (including those that warned me about a singular gradient.

I looked here but I don't see where to go form here given that my model doesn;t seem to work.

How can I model and graph this function in R?

What I tried:

LTI.func <- function(t,C,Ao,s,wd,ph) C + Ao * exp(-s*t) * cos(wd * t + ph)
mod <- nls(Y ~ LTI.func(t = I(scale(t)), C, Ao, s, wd, ph), 
           data = dat, 
           start = list(C = 1, Ao = 1, s = 1, w = 1, ph = 1))

I had no idea what starts to select, so I tried a bunch of random ones, which resulted in errors. Even when I picked starts guided by the y(t) ~ t trend I could see, I always got some kind of error:

Error in nlsModel(formula, mf, start, wts) : 
  singular gradient matrix at initial parameter estimates

Error in nls(Y ~ LTI.func(I(scale(t)), C, Ao, s, wd, ph), data = dat  : 
  singular gradient

Update:

Here is an example set of data:

dat <- data.frame(t = c(72, 25, 10, 88, 67, 63, 34, 41, 75, 13, 59, 8, 30, 52, 21),
                  Y = c(108.7, 157.5, 17.7, 175, 246.8, 233.5, 208.6, 246.5, 126.5, 
                        45.5, 214.1, 4.9, 184, 239.2, 113.3))
Z.Lin
  • 28,055
  • 6
  • 54
  • 94
theforestecologist
  • 4,667
  • 5
  • 54
  • 91
  • 1
    Which of those values are data and which are values you want to estimate? When asking for help, you should include a simple [reproducible example](https://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) with sample input and desired output that can be used to test and verify possible solutions. Show your code attempts. – MrFlick May 04 '18 at 20:57
  • Re how to plot the equation: `curve(f(x, 10, 10, 10, 10, 10))` – Frank May 04 '18 at 20:58
  • do you have the values of `y(t)` and `t`? – Onyambu May 04 '18 at 20:59
  • @Onyambu I have y(t) and t for known samples. I know nothing else (except, obviously, the graphical appearance of the y(t) ~ t relationship. – theforestecologist May 04 '18 at 21:19
  • Then you need to give some dataframe for y and t for someone to have an example of what you are talking about – Onyambu May 04 '18 at 21:21
  • 1
    My bad, guys. I'm half asleep. I've updated with code and examples – theforestecologist May 04 '18 at 21:34

1 Answers1

1

This showed up in the close queue but seemed to me to be reasonably well-formed, if poorly checked for consistency of parameter names. Here's my shot at a solution. I've tried first changing the argument to the function to x and when that didn't work tried tweaking the starting values. Eventually, I decided to scale the data argument:

LTI.func <- function(x,C,Ao,s,wd,ph) {C + Ao * exp(-s*x) * cos(wd * x + ph)}
 mod <- nls(Y ~ LTI.func(x=t , C, Ao, s, wd, ph), data=data.frame(scale(dat)), 
                     start=list(C = 0,Ao = -1,s = 1,wd = 1,ph = 0))
 mod
#-------------
Nonlinear regression model
  model: Y ~ LTI.func(x = t, C, Ao, s, wd, ph)
   data: data.frame(scale(dat))
        C        Ao         s        wd        ph 
 0.288729 -0.986426  0.517128  2.002040  2.756004 
 residual sum-of-squares: 1.53608

Number of iterations to convergence: 18 
Achieved convergence tolerance: 0.0000038776

Whether that would be a useful "solution" will require back-transforming and plotting the results against the original values or perhaps plotting the transformed coordinates against the theoretical. I wasn't surprised that the A0 value was less than zero. It certainly looked from the data that the trend was upward and exp(-s*x) would generally be downward if s*x were positive.

IRTFM
  • 258,963
  • 21
  • 364
  • 487
  • Thanks for this! I got your approach to work. However, (not due to a lack in your answer by any means) changing the parameters either does nothing to the model line (i.e., _ala_ `lines(x, predict...)`) or produces an error. I haven't quite wrapped my head around that yet. Why would changes either fail or do nothing? I threw the equation into excel to easily watch prediction values change as I changed parameters, so I know the equation itself can be more plastic. Is this an `nls` issue??? – theforestecologist May 05 '18 at 02:41
  • I don't understand what you mean by "changing the parameters". (No code.) Are we talking about the starting values? If so, why would it bother you that changing them has no effect? – IRTFM May 05 '18 at 02:51
  • Sorry, yes I meant the starting values. It only bothers me because I don;t understand why that's the case. (I'll be honest, I'm < 5 hrs deep into learning/workign with non-linear models like this, so I'm just a noob when it comes to understanding anything :p) – theforestecologist May 05 '18 at 03:06
  • 1
    The starting values are only designed to get the search started in a "region" where there might be a reasonably good possibility of finding a solution "nearby" without having to cross any sections of parameter space where the function "blows up" (goes to `Inf` or `-Inf`. So it's actually good practice to change them around and your hope is that `nls` would still end up "in the same place". – IRTFM May 05 '18 at 03:53
  • Also, I just found a neat way to scale data without having to scale the whole data.frame (which is a problem for non-numeric columns in the data): `transform(dat, t = scale(t), Y = scale(Y))`. – theforestecologist May 05 '18 at 16:46
  • For my own reference (and others) [this post](https://stackoverflow.com/q/10287545/4581200) discusses back-transforming a `scale`'d object in R. – theforestecologist May 05 '18 at 16:51
  • 1
    As long as you make an assignment that should work. The way I did it gives `nls` a scaled matrix object. Could also assign `data.frame(scale(dat))`. Also probably would be safer to use `center= FALSE` since sometimes functions will fail with negative numbers. The information about back-scaling would also be available at `?scale`. – IRTFM May 05 '18 at 17:00