[Rd] problem with function 'optimise' (PR#9438)

[Prof Brian Ripley]

That's only incidental to the story.  The problem is (code in 
src/appl/fmin.c) that at the first step symmetry (in both examples) gives 
fu=fx, and the code has

 	if (fx <= fu) {
 	    if (u < x)
 		a = u;
 	    else
 		b = u;
 	}
 	if (fu <= fx) {
...

This appears to be an error in the C translation, for the original Fortran 
does not do this.


On Thu, 4 Jan 2007, Dimitris Rizopoulos wrote:

> the issue here is a numerical one; for instance in your first example
> check:
>
> eps <- sqrt(.Machine$double.eps)
> opt1 <- optimise(ex1, lower = 0 - eps, upper = 0.5, maximum = TRUE)
>
> # or
>
> opt1 <- optimise(ex1, lower = 0 + eps, upper = 0.5, maximum = TRUE)
>
>
> which gives the correct results.
>
> I hope it helps.
>
> Best,
> Dimitris
>
> ----
> Dimitris Rizopoulos
> Ph.D. Student
> Biostatistical Centre
> School of Public Health
> Catholic University of Leuven
>
> Address: Kapucijnenvoer 35, Leuven, Belgium
> Tel: +32/(0)16/336899
> Fax: +32/(0)16/337015
> Web: http://med.kuleuven.be/biostat/
>     http://www.student.kuleuven.be/~m0390867/dimitris.htm
>
>
>
> ----- Original Message -----
> From: <karsten_krug at gmx.net>
> To: <r-devel at stat.math.ethz.ch>
> Cc: <R-bugs at biostat.ku.dk>
> Sent: Thursday, January 04, 2007 11:24 AM
> Subject: [Rd] problem with function 'optimise' (PR#9438)
>
>
>> Full_Name: Karsten Krug
>> Version: 2.4.0
>> OS: Open Suse 10.0, Windows XP
>> Submission from: (NULL) (88.134.13.50)
>>
>>
>>
>> I found a problem in the 'optimise' function for one dimensional
>> optimisation.
>>
>>
>> Example 1:
>> Try to find a maximum of the function below with the use of
>> 'optimise' in the
>> interval [0,0.5]. The function follows a parabola and has two local
>> maxima
>> located at the margins of the interval [0,0.5]. Please try the
>> following code
>> which performs this optimisation, produces a plot of the function
>> within the
>> given interval and adds the computed maximum returned by
>> 'optimise()'
>>
>> ex1 <- function(x) log(4 * x^2 - 2 * x + 0.5)
>> opt1 <- optimise(ex1, lower=0, upper=0.5, maximum=T)
>> x <- seq(0,0.5,0.001)
>> plot(x, ex1(x), type="l")
>> abline(v=opt1$maximum, col="red")
>>
>> The returned value of the maximum is 0.309017 which is definitely
>> wrong.
>>
>>
>> Example 2:
>> The next function has also two equivalent maxima this time not
>> located at the
>> margins of the interval [0, 9/41]. Thus, the maxima are real
>> analytic extrema.
>> The function 'optimise()' again returns a wrong value.
>>
>> ex2 <- function(x) log((18/41) * x - 2 * x^2) + 16 *  log(4 * x^2 -
>> (36/41) * x
>> + (9/41)) + 24 * log((23/82) +  (18/41) * x - 2 * x^2)
>> opt2 <- optimise(ex2, lower=0, upper=(9/41), maximum=T)
>> x <- seq(0,9/41,0.001)
>> plot(x,ex2(x),type="l")
>> abline(v=opt2$maximum, col="red")
>>
>>
>> The two functions are not only of theoretical interest but have to
>> be optimised
>> as likelihood functions when dealing with SNP data.
>>
>> I observed this problem on Windows XP and Open Suse 10.0. Other
>> platforms have
>> to be tested.
>>
>>
>> Some informations:
>>
>> Open Suse 10.0
>>
>> platform       i686-pc-linux-gnu
>> arch           i686
>> os             linux-gnu
>> system         i686, linux-gnu
>> status
>> major          2
>> minor          4.0
>> year           2006
>> month          10
>> day            03
>> svn rev        39566
>> language       R
>> version.string R version 2.4.0 (2006-10-03)
>>
>> Windows XP:
>>
>> platform       i386-pc-mingw32
>> arch           i386
>> os             mingw32
>> system         i386, mingw32
>> status
>> major          2
>> minor          4.0
>> year           2006
>> month          10
>> day            03
>> svn rev        39566
>> language       R
>> version.string R version 2.4.0 (2006-10-03)
>>
>>
>>
>> With best regards,
>>
>> Karsten Krug
>> Student at University of Leipzig
>>
>> ______________________________________________
>> R-devel at r-project.org mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>>
>
>
> Disclaimer: http://www.kuleuven.be/cwis/email_disclaimer.htm
>
> ______________________________________________
> R-devel at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
>

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595