[R] Robust nonlinear regression - sin(x)/x?
cstrato
cstrato at aon.at
Mon Feb 2 22:09:31 CET 2004
Dear all
Thank you very much this time for the fast response and
your many comments, and sorry for the stupid mistake.
1, The following gives the following error:
nf <- nls(z ~ a*sin(x)/(b*x), data=df,
start=list(a=0.8,b=0.9), trace = TRUE)
Error in nlsModel(formula, mf, start) : singular gradient matrix at
initial parameter estimates
2, However, as Peter Dalgaard mentioned, the following
gives the correct result:
nf <- nls(z ~ c*sin(x)/x, data=df,
start=list(c=0.5), trace = TRUE)
2.113783 : 0.5
0.3187204 : 1.022591
3, Now to the question of robust nonlinear fitting:
Let me introduce some outliers:
z1 <- z
z1[c(6,12,13,34,36,42,67,69,72,76)] <-
c(0.8,0.9,0.8,-0.5,-0.4,-0.6,0.5,0.6,0.8,0.7)
plot(x,z1)
df1 <- as.data.frame(cbind(x,z1))
Now, the fit gives:
nf1 <- nls(z1 ~ c*sin(x)/x, data=df1,
start=list(c=0.5), trace = TRUE)
4.814774 : 0.5
2.072135 : 1.145962
The true result should be c=1.0
fitting w/o outliers gives c=1.023
fitting with outliers gives c=1.146
Can this fit considered to be robust?
Best regards
Christian
Peter Dalgaard wrote:
> cstrato <cstrato at aon.at> writes:
>
>
>>Dear all
>>
>>Since I did not receive any answer to my general question (?),
>>let me ask a concrete question:
>>
>>How can I fit the simple function y = a*sin(x)/b*x?
>>
>>This is the code that I tried, but nls gives an error:
>>
>>x <- seq(1,10,0.1)
>>y <- sin(x)/x
>>plot(x,y)
>>z <- jitter(y,amount=0.1)
>>plot(x,z)
>>df <- as.data.frame(cbind(x,z))
>>nf <- nls(z ~ a*sin(x)/b*x, data=df,
>> start=list(a=0.8,b=0.9), trace = TRUE)
>>
>>I have followed the Puromycin sample which works fine:
>>Pur.wt <- nls(rate ~ (Vm * conc)/(K + conc), data = Treated,
>> start = list(Vm = 200, K = 0.1), trace = TRUE)
>>
>>Do I make some mistake or is it not possible to fit sin(x)/x?
>
>
> The expression only depends on a/b so you cannot estimate both.
>
> Besides, you need to check up on operator precedence and
> associativity: What you wrote is equivalent to a*sin(x)*x/b.
>
More information about the R-help
mailing list