Thanks a lot Ravi for the clarification and it makes sense - bug in my understanding acknowledged.
I'll change the objective function as suggested.
Thanks,
Rainer
Ravi Varadhan <ravi.varad...@jhu.edu> writes:
> Dear Rainer,
> This is NOT a bug in auglag. I already mentioned that auglag() can
> work with infeasible starting values, which also implies that the
> function must be evaluable at infeasible values. A simple solution to
> your problem would be to fix up your objective function such that it
> evaluates to `Inf' or some large value, when the parameter values are
> not in the constrained domain. constrOptim.nl() is a barrier method
> so it forces the initial value and the subsequent iterates to be
> feasible.
> Best,
> Ravi
> ________________________________________
> From: Rainer M Krug <rai...@krugs.de>
> Sent: Tuesday, October 6, 2015 9:20 AM
> To: Ravi Varadhan
> Cc: 'r-help@r-project.org'
> Subject: Bug in auglag?
>
> Hi Ravi,
>
> I would like come back to your offer. I have a problem which possibly is
> caused by a bug or by something I don't understand:
>
> My function to be minimised is executed even when an element in hin() is
> negative.
>
> My hin looks as follow:
>
> hinMahat <- function(x, hauteur, na, zjoint, y, LAI, ...) {
> if (x[1] < 0) {
> cat(names(list(...)), "\n")
> cat(..., "\n")
> cat(x, "|", hauteur, LAI, y, "\n")
> }
>
> h <- rep(NA, 8)
> if (!missing(na)) {
> x <- c(na, x )
> }
> if (!missing(y)) {
> x <- c(x, y)
> }
> if (!missing(zjoint)) {
> x <- c(x[1], zjoint, x[2])
> }
>
> ##
> dep <- hauteur * (0.05 + LAI^0.02 / 2) + (x[3] - 1)/20
> h[1] <- dep
> h[2] <- hauteur - dep
> ## if (h[2]==0) {
> ## h[2] <- -1
> ## }
> ##
> z0 <- hauteur * (0.23 + LAI^0.25 / 10) + (x[3] - 1)/67
> h[3] <- z0
> ## if (h[3]==0) {
> ## h[3] <- -1
> ## }
> h[4] <- hauteur - z0
> ##
> h[5] <- x[1]
> ##
> h[6] <- x[2]
> h[7] <- hauteur - x[2]
> ##
> h[8] <- hauteur - dep - z0
> if (any(h<=0)) {
> cat(h, "\n")
> cat("\n")
> }
> return(h)
> }
>
> the x contains up to three elements: c(na=, zjoint=, y=) and I fit these
> three, unless one or two are specified explicitely.
>
> The values going into hin are:
>
> ,----
> | ... (z u ua za z0sol )
> | 3 11 17 23 29 37 0.315 0.422 0.458 0.556 1.567 1.747 1.747 37 0.001
> |
> | x(na, zjoint): -8.875735 24.51316
> | hauteur: 28
> | na: 8.1
> | y: 3
> |
> | the resulting hin() is:
> | 16.09815 11.90185 11.19352 16.80648 -8.875735 24.51316 3.486843 0.708335
> `----
>
>
> Which is negative in element 5 as x[2]=na is negative.
>
> So I would expect that the function fn is not evaluated. But it is, and
> raises an error:
>
> ,----
> | Error in wpLELMahat(z = z, ua = ua, na = ifelse(missing(na), par[1], na), :
> | na has to be larger or equal than zero!
> `----
>
> Is this a misunderstanding on my part, or is it an error in the function
> auglag?
>
>
> Below is the function which is doing the minimisation.
>
> If I replace auglag() with constrOptim.nl(), the optimisation is working
> as expected.
>
> So I think this is a bug in auglag?
>
> Let me know if you need further information.
>
> Cheers,
>
> Rainer
>
> --8<---------------cut here---------------start------------->8---
> fitAuglag.wpLEL.mahat.single <- function(
> z,
> u,
> LAI,
> initial = c(na=9, zjoint=0.2*2, y=3),
> na, zjoint, y,
> h = 28,
> za = 37,
> z0sol = 0.001,
> hin,
> ...
> ) {
> if (missing(hin)) {
> hin <- hinMahat
> }
>
> wpLELMin <- function(par, na, zjoint, y, z, u, ua, hauteur, za, z0sol,
> LAI) {
> result <- NA
> try({
> p <- wpLELMahat(
> z = z,
> ua = ua,
> na = ifelse(missing(na), par[1], na),
> zjoint = ifelse(missing(zjoint), par[2], zjoint),
> h = hauteur,
> za = za,
> z0sol = z0sol,
> LAI = LAI,
> y = ifelse(missing(y), par[3], y)
> )
> result <- sum( ( (p$u - u)^2 ) / length(u) )
> },
> silent = FALSE
> )
> ## cat("From wpLELMin", par, "\n")
> return( result )
> }
>
> ua <- u[length(u)]
> result <- list()
> result$method <- "fitAuglag.wpLEL.mahat.single"
> result$initial <- initial
> result$dot <- list(...)
> result$z <- z
> result$u <- u
>
> result$fit <- auglag(
> par = initial,
> fn = wpLELMin,
> hin = hin,
> na = na,
> zjoint = zjoint,
> y = y,
> ##
> z = z,
> u = u,
> ua = ua,
> hauteur = h,
> za = za,
> z0sol = z0sol,
> LAI = LAI,
> ...
> )
> result$wp <- wpLELMahat(
> z = z,
> ua = ua,
> na = ifelse ( missing(na), result$fit$par["na"], na),
> zjoint = ifelse ( missing(zjoint), result$fit$par["zjoint"], zjoint),
> h = h,
> za = za,
> z0sol = z0sol,
> LAI = LAI,
> y = ifelse ( missing(y), result$fit$par["y"], y)
> )
>
> class(result) <- c(class(result), "wpLELFit")
> return(result)
> }
> #+end_src--8<---------------cut here---------------end--------------->8---
>
>
>
> Ravi Varadhan <ravi.varad...@jhu.edu> writes:
>
>> I would recommend that you use auglag() rather than constrOptim.nl()
>> in the package "alabama." It is a better algorithm, and it does not
>> require feasible starting values.
>> Best,
>> Ravi
>>
>> -----Original Message-----
>> From: Rainer M Krug [mailto:rai...@krugs.de]
>> Sent: Thursday, October 01, 2015 3:37 AM
>> To: Ravi Varadhan <ravi.varad...@jhu.edu>
>> Cc: 'r-help@r-project.org' <r-help@r-project.org>
>> Subject: Re: optimizing with non-linear constraints
>>
>> Ravi Varadhan <ravi.varad...@jhu.edu> writes:
>>
>>> Hi Rainer,
>>> It is very simple to specify the constraints (linear or nonlinear) in
>>> "alabama" . They are specified in a function called `hin', where the
>>> constraints are written such that they are positive.
>>
>> OK - I somehow missed the part that, when the values x are valid,
>>> i.e. in the range as defined by the conditions, the result of hin(x)
>>> that they are all positive.
>>
>>> Your two nonlinear constraints would be written as follows:
>>>
>>> hin <- function(x, LAI) {
>>> h <- rep(NA, 2)
>>> h[1] <- LAI^x[2] / x[3] + x[1]
>>> h[2] <- 1 - x[1] - LAI^x[2] / x[3]
>>> h
>>> }
>>
>> Makes perfect sense.
>>
>>>
>>> Please take a look at the help page. If it is still not clear, you can
>>> contact me offline.
>>
>> Yup - I did. But I somehow missed the fact stated above.
>>
>> I am using constrOptim() and constrOptim.nl() for a paper and am
>>> compiling a separate document which explains how to get the
>>> constraints for the two functions step by step - I will make it
>>> available as a blog post and a pdf.
>>
>> I might have further questions concerning the different fitting
>>> functions and which ones are the most appropriate in my case.
>>
>> Thanks a lot,
>>
>> Rainer
>>
>>
>>> Best,
>>> Ravi
>>>
>>> Ravi Varadhan, Ph.D. (Biostatistics), Ph.D. (Environmental Engg)
>>> Associate Professor, Department of Oncology Division of Biostatistics
>>> & Bionformatics Sidney Kimmel Comprehensive Cancer Center Johns
>>> Hopkins University
>>> 550 N. Broadway, Suite 1111-E
>>> Baltimore, MD 21205
>>> 410-502-2619
>>>
>>>
>>> [[alternative HTML version deleted]]
>>>
>>
>> --
>> Rainer M. Krug
>> email: Rainer<at>krugs<dot>de
>> PGP: 0x0F52F982
>>
>
> --
> Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation Biology,
> UCT), Dipl. Phys. (Germany)
>
> Centre of Excellence for Invasion Biology
> Stellenbosch University
> South Africa
>
> Tel : +33 - (0)9 53 10 27 44
> Cell: +33 - (0)6 85 62 59 98
> Fax : +33 - (0)9 58 10 27 44
>
> Fax (D): +49 - (0)3 21 21 25 22 44
>
> email: rai...@krugs.de
>
> Skype: RMkrug
>
> PGP: 0x0F52F982
>
--
Rainer M. Krug
email: Rainer<at>krugs<dot>de
PGP: 0x0F52F982
signature.asc
Description: PGP signature
______________________________________________ R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see 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.
