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.0225913, 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.145962The 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 <[EMAIL PROTECTED]> 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.
______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
