Skip to content

fitting "power model" in nls()

5 messages · Milton Cezar Ribeiro, Rolf Turner, Joerg van den Hoff +1 more

Rolf Turner
# 
Rule number 1:  Read the help for the function you are using.
You must supply starting values for the fit --- which the code you gave
doesn't do.

Rule number 2:  Don't use nls()!  Endless grief results.

Instead try:

foo <- function(par,x,y){
  Const <- par[1]
  B <- par[2]
  A <- par[3]
  sum((y-(Const+B*(x^A)))^2)
  }

fit <- optim(c 
(1,2.5,0.2),foo,x=area,y=richness,method="BFGS",control=list 
(maxit=1000))

plot(area,richness)
ccc <- fit$par
curve(ccc[1]+ccc[2]*(x^ccc[3]),from=0,to=2000,add=TRUE,col="red") #  
Fit doesn't look too bad.

HTH.

		cheers,

			Rolf Turner
On 3/12/2007, at 8:08 AM, Milton Cezar Ribeiro wrote:

            

      
      
      
    

        

      
      
    

        
######################################################################
Attention: 
This e-mail message is privileged and confidential. If you are not the 
intended recipient please delete the message and notify the sender. 
Any views or opinions presented are solely those of the author.

This e-mail has been scanned and cleared by MailMarshal 
www.marshalsoftware.com
######################################################################

  
  


  


      

    

    

      
      

        


  

        

      Joerg van den Hoff
    

    

      
      #
    

  


  

      
On Sun, Dec 02, 2007 at 11:08:01AM -0800, Milton Cezar Ribeiro wrote:

            

      
      
      
    

        
for easier notation, let y=richness, x=area, C=const in the following.

then 

nls(y~B*x^A + C, start = c(A=3.2, B=0.002, C=0))

converges alright. where's the problem (apart from this being not a very good
model for the data)? the critical point is to provide some reasonable estimate
of the parameters as starting values.

to get reasonable start values, I'd use:


y = B*x^A + C --> log(y-C) = log(B) + A*log(x) --> ly = b + a*lx,

estimate C from the x -> 0 asymptotic value (approx. 0)

and use lm(ly~lx) which yields a and b estimates which you could use in `nls'.

and, contrary to other assessments you've received, I definitely would prefer `nls'
for least squares fitting instead of using `optim' or other general minimization routines.

hth

joerg

  
  


  


      

    

    

      
      

        


  

        

      Gabor Grothendieck
    

    

      
      #
    

  


  

      
Is that really the model we want?  When we have problems sometimes
its just a sign that the model is not very good in the first place.

plot(richness ~ area)

shows most of the points crowded the left and just a few points out to
the right.  This
does not seem like a very good pattern for model fitting.

plot(richness ~ log(area))
plot(log(richness) ~ log(area))

both look nicer.

        
On Dec 2, 2007 2:08 PM, Milton Cezar Ribeiro <milton_ruser at yahoo.com.br> wrote:

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Rolf Turner
    

    

      
      #
    

  


  

      
On 3/12/2007, at 9:26 AM, Joerg van den Hoff wrote:

        
<snip>

            

      
      
      
    

        
Clearly you are far cleverer than I at getting nls() to work.

I attempted to replicate your recipe ***exactly***, right down to the  
very last
detail of notation and the ordering of terms in the expression:

  > x <- area
  > y <- richness
  > ml2 <- nls(y ~ B*x^A + C,start=c(A=3.2,B=0.002,C=0))

and got the following error:

Error in nls(y ~ B * x^A + C, start = c(A = 3.2, B = 0.002, C = 0)) :
   singular gradient

	cheers,

		Rolf Turner

P.S.:

Version information:

 > version
                _
platform       i386-apple-darwin8.10.1
arch           i386
os             darwin8.10.1
system         i386, darwin8.10.1
status
major          2
minor          6.0
year           2007
month          10
day            03
svn rev        43063
language       R
version.string R version 2.6.0 (2007-10-03)

######################################################################
Attention:\ This e-mail message is privileged and confid...{{dropped:9}}