On May 4, 2010, at 5:25 PM, ivo welch wrote:
thank you, david and gabor. very much appreciated. I should have
thought of setting the seed. this was only an example, of course.
alas, such intermittent errors could still be of concern to me,
because I need to simulate this nls() to find out its properties under
the NULL, so I can't easily tolerate errors. fortunately, I had the
window still open, so getting my y's out was easy, and the rounded
figures produce the same nls error.
There would of course be try()
But ... Again, No error on my device:
> xy <- scan()
1: 1 5.017
3: 2 7.993
5: 3 11.014
7: 4 13.998
9: 5 17.003
11: 6 19.977
13: 7 23.011
15: 8 25.991
17: 9 29.003
19: 10 32.014
21: 11 31.995
23: 12 32.004
25: 13 32.012
27: 14 31.994
29: 15 31.998
31: 16 32.000
33: 17 32.009
35: 18 31.995
37: 19 32.000
39: 20 31.982
41:
Read 40 items
> xym <-matrix(xy, ncol=2, byrow=TRUE)
> colnames(xym)<-c("x","y")
> r1= nls( y~ a+b*pmin(c,x), start=list(a=2, b=3, c=10), trace=TRUE )
0.001505770 : 2 3 10
0.001361737 : 2.003063 2.999924 10.000135
Plotting predict(r1) confirmed a very close fit.
> sessionInfo()
R version 2.10.1 RC (2009-12-09 r50695)
x86_64-apple-darwin9.8.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/C/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] misc3d_0.7-0 rgl_0.91 spatstat_1.18-3 deldir_0.0-12
[5] mgcv_1.6-1 ks_1.6.12 mvtnorm_0.9-9
KernSmooth_2.23-3
[9] lattice_0.18-3
loaded via a namespace (and not attached):
[1] grid_2.10.1 Matrix_0.999375-38 nlme_3.1-96
tools_2.10.1
cbind(x,round(y,3))
x y
[1,] 1 5.017
[2,] 2 7.993
[3,] 3 11.014
[4,] 4 13.998
[5,] 5 17.003
[6,] 6 19.977
[7,] 7 23.011
[8,] 8 25.991
[9,] 9 29.003
[10,] 10 32.014
[11,] 11 31.995
[12,] 12 32.004
[13,] 13 32.012
[14,] 14 31.994
[15,] 15 31.998
[16,] 16 32.000
[17,] 17 32.009
[18,] 18 31.995
[19,] 19 32.000
[20,] 20 31.982
r1= nls( y~ a+b*pmin(c,x), start=list(a=2, b=3, c=10), trace=TRUE )
0.002138 : 2 3 10
0.002117 : 2.004 3.000 9.999
0.002113 : 2.006 2.999 10.001
0.002082 : 2.005 2.999 10.000
0.002077 : 2.005 2.999 10.000
0.002077 : 2.005 2.999 10.000
Error in nls(y ~ a + b * pmin(c, x), start = list(a = 2, b = 3, c =
10), :
step factor 0.000488281 reduced below 'minFactor' of 0.000976562
I really don't care about this example, of course---only about
learning how to avoid nls() from dying on me. so, any advice would be
appreciated.
regards,
/iaw
----
Ivo Welch (ivo.we...@brown.edu, ivo.we...@gmail.com)
On Tue, May 4, 2010 at 3:59 PM, David Winsemius <dwinsem...@comcast.net
> wrote:
On May 4, 2010, at 3:52 PM, ivo welch wrote:
thank you, david. indeed. works great (almost). an example for
anyone else googling this in the future:
x=1:20
y= 2+3*ifelse(x>10, 10, x)+rnorm(20,0,0.01)
r1= nls( y~ a+b*pmin(c,x), start=list(a=2, b=3, c=10), trace=TRUE )
0.002142 : 2 3 10
0.002115 : 2.004 3.000 10.000
0.002114 : 2.006 2.999 10.001
0.002084 : 2.005 2.999 10.000
...
0.002079 : 2.005 2.999 10.000
Error in nls(y ~ a + b * pmin(c, x), start = list(a = 2, b = 3, c
= 10),
:
step factor 0.000488281 reduced below 'minFactor' of 0.000976562
strange error, but unrelated to my question. will figure this one
out
next.
I get no error. May be difficult to sort out unless you can
reproduce after
setting a random seed.
x=1:20
y= 2+3*ifelse(x>10, 10, x)+rnorm(20,0,0.01)
r1= nls( y~ a+b*pmin(c,x), start=list(a=2, b=3, c=10), trace=TRUE )
0.001560045 : 2 3 10
0.001161253 : 2.003824 2.998973 10.000388
0.001161253 : 2.003824 2.998973 10.000388
--
David.
regards,
/iaw
On Tue, May 4, 2010 at 3:40 PM, David Winsemius <dwinsem...@comcast.net
>
wrote:
On May 4, 2010, at 3:33 PM, ivo welch wrote:
Dear R experts---I would like to estimate a non-linear least
squares
expression that looks something like
y ~ a+b*min(c,x)
where a, b, and c are the three parameters. how do I define a min
function in the formula language of R? advice appreciated.
?pmin
sincerely,
/iaw
David Winsemius, MD
West Hartford, CT
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
