On 5/17/2006 11:07 AM, Martin Maechler wrote:
>>>>>> "Duncan" == Duncan Murdoch <[EMAIL PROTECTED]>
>>>>>> on Tue, 16 May 2006 08:34:06 -0400 writes:
>
> Duncan> On 5/16/2006 4:56 AM, [EMAIL PROTECTED]
> Duncan> wrote:
> >> Probably I included too much at once in my bug report. I
> >> can live with an unfulfilled wishlist and thank you for
> >> thinking about it. The "badly-behaved" function is just
> >> an example to demonstrate the bug I reported. I think it
> >> is a bug if optim returns (without any warning) an
> >> unmatching pair of par and value: f(par) != value. And it
> >> is easily fixed.
>
> >> Andreas
>
> Duncan> I agree with you that on return f(par) should be
> Duncan> value. I agree with Brian that changes to the
> Duncan> underlying strategy need much more thought.
>
> I agree (to both).
> However, isn't Andreas' patch just fixing the problem
> and not changing the underlying strategy at all?
> [No, I did not study the code in very much detail ...]
Brian and I only quoted part of his message. The patch we quoted isn't
bad, but I'm not sure it's the best: in particular, with the patch
optim() returns a function value that is larger than f at the starting
value (see below). I think this means it would be better to change
optim$par rather than changing optim$value to achieve consistency, but a
quick look at optim.c made me think it would take more time than I had
to do this without messing up something else.
I don't think I'll have a chance to look at this before 2.3.1, so if
nobody else takes it on, I'd prefer to leave this as an unresolved bug
report for now.
>
> Martin Maechler
>
> >> Prof Brian Ripley wrote:
> >>
> >>> [Sorry for the belated reply: this came in just as I was leaving for a
> >>> trip.]
> >>>
> >>> I've checked the original source, and the C code in optim does
> >>> accurately reflect the published algorithm.
> >>>
> >>> Since your example is a discontinuous function, I don't see why you
> >>> expect CG to work on it. John Nash reports on his extensive
> >>> experience that method 3 is the worst, and I don't think we should let
> >>> a single 2D example of a badly-behaved function override that.
> >>>
> >>> Note that no other optim method copes with the discontiuity here: had
> >>> your reported that it would have been clear that the problem was with
> >>> the example.
> >>>
> >>> On Fri, 21 Apr 2006, [EMAIL PROTECTED] wrote:
> >>>
> >>>> Dear R team,
> >>>>
> >>>> when using optim with method "CG" I got the wrong $value for the
> >>>> reported $par.
> >>>>
> >>>> Example:
> >>>> f<-function(p) {
> >>>> if (!all(p>-.7)) return(2)
> >>>> if (!all(p<.7)) return(2)
> >>>> sin((p[1])^2)*sin(p[2])
> >>>> }
> >>>> optim(c(0.1,-0.1),f,method="CG",control=list(trace=0,type=1))
> >>>> $par 19280.68 -10622.32
> >>>> $value -0.2346207 # should be 2!
> >>>> optim(c(0.1,-0.1),f,method="CG",control=list(trace=0,type=2))
> >>>> $par 3834.021 -2718.958
> >>>> $value -0.0009983175 # should be 2!
I think this is f(0.1, -0.1), so really $par should be 0.1, -0.1 in this
case. In the one above, it appears to have made a little progress
before it went off track, but -0.234 is better than 2, so it should be
returned if it's really an f(p) value.
Duncan Murdoch
> >>>>
> >>>> Fix:
> >>>> --- optim.c (Revision 37878)
> >>>> +++ optim.c (Arbeitskopie)
> >>>> @@ -970,7 +970,8 @@
> >>>> if (!accpoint) {
> >>>> steplength *= stepredn;
> >>>> if (trace) Rprintf("*");
> >>>> - }
> >>>> + } else
> >>>> + *Fmin = f;
> >>>> }
> >>>> } while (!(count == n || accpoint));
> >>>> if (count < n) {
> >>>>
> >>>> After fix:
> >>>> optim(c(0.1,-0.1),f,method="CG",control=list(trace=0,type=1))
> >>>> $par 0.6993467 -0.4900145
> >>>> $value -0.2211150
> >>>> optim(c(0.1,-0.1),f,method="CG",control=list(trace=0,type=2))
> >>>> $par 3834.021 -2718.958
> >>>> $value 2
> >>>>
> >>>> Wishlist:
> >>>
> >> [wishlist deleted]
> >>
> >>
>
> Duncan> ______________________________________________
> Duncan> R-devel@r-project.org mailing list
> Duncan> https://stat.ethz.ch/mailman/listinfo/r-devel
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
