[R] Robust nonlinear regression - sin(x)/x?

[cstrato]

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.
>