[R] Newbie Question: Shifted Power Fit?

[Douglas Bates]

Bill Simpson <wsi at gcal.ac.uk> writes:

> > I have multiple sets of 2D data that are coming from
> > distributions of the form:
> > 
> > y = A(x-C)^B     (Eq.1)
> > 
> > I am trying to estimate for each set the best values 
> > of A, B, and C so that Eq.1 will be the best fit for
> > the data. I guess that it should be easy to do, but
> > I lack the experience :(
> Do this:
>  library(nls)
>  ?nls
> Here is an example
> x<-3:13
> y<-2*(x-3)^1.5  #fake perfect data
> fit<-nls(y~a*(x-b)^c,start=list(a=2,b=3,c=1.5))
> 
> This gives:
> Error in numericDeriv(form[[3]], names(ind), env) : 
>         Missing value or an Infinity produced when evaluating the model
> Maybe others here can say what is wrong.

Trying to evaluate the exponential of a negative number?  Think of how
the expression (x-b)^c for non-integer c is evaluated when x = 3 and b
= 3 - eps. 

By the way, it is a very bad idea to use c as an identifier in R.

> I personally tend to use nlm() and minimize the sum of squared errors...
> 
> Once it works do summary(fit) to see the fit results.

A more realistic example would include noise in the observations and
not have the simulated value of b equal to (or greater than) one of
the x values. 

You can use the "plinear" algorithm in nls for this example because
the parameter a appears linearly in the model.

> x <- 4:13
> y <- 2*(x-3)^1.5 + rnorm(length(x), 0, 0.1)
> fm1 <- nls(y ~ (x - b)^pow, start = c(b = 2.5, pow = 1.0),
+    algorithm = "plinear", trace = TRUE)
362.419 : 2.500000 1.000000 5.153966 
20.87003 : 1.695023 1.671158 1.080315 
1.139375 : 3.251409 1.466385 2.252482 
0.07488055 : 3.018406 1.496594 2.022821 
0.07381498 : 3.012420 1.497047 2.018591 
0.07381497 : 3.012403 1.497049 2.018575 
> summary(fm1)

Formula: y ~ (x - b)^pow

Parameters:
     Estimate Std. Error t value Pr(>|t|)    
b     3.01240    0.04386   68.69 3.64e-11 ***
pow   1.49705    0.01180  126.88 4.98e-13 ***
.lin  2.01858    0.06621   30.49 1.05e-08 ***
---
Signif. codes:  0 `***' 0.001 `**' 0.01 `*' 0.05 `.' 0.1 ` ' 1 

Residual standard error: 0.1027 on 7 degrees of freedom

Correlation of Parameter Estimates:
           b     pow
pow  -0.9482        
.lin  0.9704 -0.9963
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._