[R] R formula language---a min and max function?
Gabor Grothendieck
ggrothendieck at gmail.com
Wed May 5 01:39:41 CEST 2010
For fixed c, y is linear in pmin(c, x) so the first two statements
find c and the next solves the remaining linear problem
f <- function(c) sum((y - fitted(lm(y ~ pmin(c, x))))^2)
fit.c <- optimize(f, c(0, max(DF$x))); fit.c
lm(y ~ pmin(fit.c$minimum, x))
as a function of c
On Tue, May 4, 2010 at 5:25 PM, ivo welch <ivo.welch at gmail.com> 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.
>
>> 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.welch at brown.edu, ivo.welch at gmail.com)
>
>
>
> On Tue, May 4, 2010 at 3:59 PM, David Winsemius <dwinsemius at 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 <dwinsemius at 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
>>>>>
>>>>> ______________________________________________
>>>>> R-help at 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.
>>>>
>>>> David Winsemius, MD
>>>> West Hartford, CT
>>>>
>>>>
>>
>> David Winsemius, MD
>> West Hartford, CT
>>
>>
>
> ______________________________________________
> R-help at 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.
>
More information about the R-help
mailing list