[R] manipulating text to generate different formulas to use in nls()

Héctor Villalobos hvillalo at ipn.mx
Mon Aug 10 22:54:01 CEST 2009


Hello,

In doing a series of non-linear estimations of a function which is a sum of a varying number 
of sinusoids, I would like to "autogenerate" the arguments needed by nls() depending on that 
number.

For example, when there are two sinusoids: 

  > nls( y ~ mu + A1 * cos(2*pi*f1*x - P1) + A2 * cos(2*pi*f2*x - P2), data = some.xy.data, 
            start = list( mu=some.value0, A1=some.value1, P1=some.value2, A2=some.value3,
                         P2=some.value4, f1=some.value5, f2=some.value6  )  )

 Adding a third sinusoid, the model and the starting parameters needed become:

  > nls( y ~ mu + A1 * cos(2*pi*f1*x - P1) + A2 * cos(2*pi*f2*x - P2) + A3 * cos(2*pi*f3*x - P3), 
            data = some.xy.data, start = list( mu=some.value0, A1=some.value1,
            P1=some.value2, A2=some.value3, P2=some.value4, A3=some.value5,
            P3=some.value6, f1=some.value7, f2=some.value8, f3=some.value9 )


So far, I managed to create the "building blocks" I need but when it comes to defining their 
number I still haven't found a solution (see the following code).

##  some fake starting values
par.ini <- data.frame(period=1:5, invR=1:5, a0=1:5, Amplitude=1:5, Phase=1:5)

  n <- nrow(par.ini) - 1                        # I'm not interested in the last harmonic
  inip <- par.ini[1:n, c(1, 3:5)]             # neither in the second column
    amps <- paste("A", 1:n, sep="")
    phas <- paste("P", 1:n, sep="")
    freqs <- paste("f", 1:n, sep="")

 ( premodl <- paste(amps, " * cos(2*pi*", freqs, "*x - ", phas, ")", sep="") )


## In the two following instructions, "modl" and "initP" are defined by hand. This is what I 
want to generate "automagically" depending on the number of harmonics (n).

 ( modl <- paste("y ~ mu", premodl[1],premodl[2], premodl[3], premodl[4], sep=" + ") ) 

 ( initP <- list(mu=sum(inip$a0), 
          A1=inip$Amplitude[1], P1=inip$Phase[1], 
          A2=inip$Amplitude[2], P2=inip$Phase[2], 
          A3=inip$Amplitude[3], P3=inip$Phase[3],           
          A4=inip$Amplitude[4], P4=inip$Phase[4], 
          f1=1/inip$period[1], f2=1/inip$period[2],
          f3=1/inip$period[3], f4=1/inip$period[4])  ) 

## not run. These two will be used in
  res <- nls(modl, data=data.frame(y, x), start=initP)
 


Thank you in advance for any insigth.

Héctor


 
-- 
Héctor Villalobos <hvillalo at ipn.mx> 
 CICIMAR - IPN
 A.P. 592. Col. Centro 
 La Paz, Baja California Sur, MÉXICO. 23000
 Tels. (+52 612) 122 53 44; 123 46 58; 123 47 34  ext. 82425
 Fax.  (+52 612) 122 53 22




More information about the R-help mailing list