1

I am trying to run the following code but it gives me:

Error in qr.default(.swts * attr(rhs, "gradient")) : NA/NaN/Inf in foreign function call (arg 1) In addition: Warning message: In log(.expr4) : NaNs produced

Can you please help me about? Thanks!

model <- deriv( ~ c*(1+b*(q-1)*t)^(1/(1-q)),  c("c", "b", "q"), function (t, c, b,q){})

nls(Frequency ~ model(t, c, b, q), data=DF,start=list(c = 1, b = 1.5, q =0.5))

Following you can see a part of data, for which I am trying to fit a q-exponential distribution functin as I explained above. I am uiung nls function in R to obtain estimates (q-exponential) for the given data.

t Frequency 0 195746 1 93938 2 53181 3 31853 4 19856 5 12182 6 7847 7 5459 8 4325 9 3203 10 2750

13554N
  • 37
  • 9
  • 1
    Welcome to Stack Overflow! Can you please include data that will provide us with a [reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) ? – Ben Bolker Dec 19 '16 at 01:18
  • 1
    just guessing, but for your starting parameters, the base of the exponent is negative for t>4/3; any computation that involves fractional powers of negative values will give `NaN` ... – Ben Bolker Dec 19 '16 at 02:56
  • can you please edit your question to include the data (rather than having it just in a comment)? – Ben Bolker Dec 19 '16 at 16:38

1 Answers1

1

tl;dr you need to work harder to think about/find reasonable starting values.

Set up data and model:

 dd <- data.frame(t=0:10,
            Frequency=c(195746,93938,53181,31853,19856,12182,
                        7847,5459,4325,3203,2750))
model <- deriv( ~ c*(1+b*(q-1)*t)^(1/(1-q)),
               c("c", "b", "q"), function (t, c, b,q){})

If you just evaluate model() at the initial parameters you can see that some of the derivatives with respect to q are NaN (it is always worthwhile to try evaluating the objective function and gradient at the starting values to make sure they make sense). In addition to this, if you look at the values of the objective function you can see that they're nowhere close to the data (they increase from 1 to 42, whereas the data decrease from 200000 to 3000 ...)

par0 <- list(c = 1, b = 1.5, q =0.5)
with(par0,model(dd$t,c,b,q))
##  [1]  1.0000  0.0625  0.2500  1.5625  4.0000  7.5625 12.2500 18.0625 25.0000
## [10] 33.0625 42.2500
## attr(,"gradient")
##             c     b         q
##  [1,]  1.0000  0.00 0.0000000
##  [2,]  0.0625 -0.25 0.4034264
##  [3,]  0.2500  1.00       NaN
##  [4,]  1.5625  3.75       NaN
...

The problem is that raising a negative value to a fractional power gives NaN. When q is <1, the base will go negative when t is sufficiently large ...

with(par1,model(dd$t,c,b,q))

I tried setting q=0 (so that the base value of the exponent is 1), but that still encounters the same problem.

What if we start from values where q-1 is >0 so that the base will be positive?

par2 <- list(c = 1, b = 1.5, q =2)
nls(Frequency ~ model(t, c, b, q), data=dd,start=par2)

Error in nls(Frequency ~ model(t, c, b, q), data = dd, start = par2) : step factor 0.000488281 reduced below 'minFactor' of 0.000976562

OK, that's getting a little better but still failing. What if we try to start on a more reasonable scale, i.e. set c to 100000?

par3 <- list(c = 1e5, b = 1.5, q =2)
fit3 <- nls(Frequency ~ model(t, c, b, q), data=dd,start=par3)

Plot results:

par(las=1,bty="l")
plot(Frequency/1000~t,data=dd)
newdat <- data.frame(t=seq(0,10,length=51))
lines(newdat$t,predict(fit3,newdata=newdat)/1000)

enter image description here

Ben Bolker
  • 211,554
  • 25
  • 370
  • 453
  • As I told I am fitting a q-exponential function to my data, and for this 00. I am trying to figure out your answer. In case I will write you back. Thank you very much. – 13554N Dec 19 '16 at 17:21
  • It seem it is working well. Another question is that how can I have a self-starting estimate, I mean with of providing staring values for the parameters? Thank you! – 13554N Dec 19 '16 at 17:23
  • 1
    In order to come up with a self-starting estimate you'd need to come up with your own algorithm for finding approximate initial starting values (e.g. see the existing self-start functions, `apropos("^SS[[:lower:]]")`). You could ask another question along those lines (but you should see how far you can get on your own first). If this answer was useful you can upvote it (if you have sufficient reputation), and in any case if it answers your question satisfactorily you are encouraged to click the check-mark to accept it. – Ben Bolker Dec 19 '16 at 17:25
  • Well, I developed two other nonlinear models. Now, I want to compare three models fitness and see which one is a better fit. Does the "residual sum-of-squares" is a good index in this regard, (the lower is the better fit)? Or is there some technique to do these kind of comparisons for nonlinear regression models? Thanks! – 13554N Dec 20 '16 at 00:00
  • residual sum of squares is reasonable if the models have the same number of parameters, otherwise it gets complicated. That would definitely be a [CrossValidated](http://stats.stackexchange.com) question rather than a SO question though ... – Ben Bolker Dec 20 '16 at 14:49
  • Regarding starting values, since the objective function is linear in `c` and `c*b` you could use `algorithm = "plinear"` in which case you only need starting values for the nonlinear parameters, viz. `q`. Using a grid of possible `q` values choose the best one. You could use `nls2` in the nls2 package to easily do this. – G. Grothendieck Dec 21 '16 at 13:32
  • @BenBolker I asked my question in CrossValidated and they raised the question that my models are not nested therefore I cannot use anova (), however I could not figure out how can I compare the model is built above with following model: model3 <- deriv( ~ exp(-beta1*t)/(beta2+beta3*t),c("beta1", "beta2", "beta3"), function (t, beta1, beta2,beta3){}) par3 <- list(beta1 = .0100, beta2 = .0002, beta3 =0.0002) fit3 <- nls(Frequency ~ model3(t, beta1,beta2,beta3), data=dd,start=par3). My question is that aren't these models nested? and how can I compare them? – 13554N Dec 21 '16 at 23:45
  • @G.Grothendieck, thanks a lot, I will try the method you mentioned. – 13554N Dec 21 '16 at 23:47