[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