Message-ID: <401EBC8B.3010908@aon.at>
Date: 2004-02-02T21:09:31Z
From: cstrato
Subject: Robust nonlinear regression - sin(x)/x?
In-Reply-To: <x2isip45lv.fsf@biostat.ku.dk>
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.
>