[R] Bug in auglag?

Rainer M Krug Rainer at krugs.de
Tue Oct 6 15:27:58 CEST 2015


Please ignore - list members - accidentally CCd.

Rainer


Rainer M Krug <Rainer at krugs.de> writes:

> 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.varadhan at 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:Rainer at krugs.de] 
>> Sent: Thursday, October 01, 2015 3:37 AM
>> To: Ravi Varadhan <ravi.varadhan at jhu.edu>
>> Cc: 'r-help at r-project.org' <r-help at r-project.org>
>> Subject: Re: optimizing with non-linear constraints
>>
>> Ravi Varadhan <ravi.varadhan at 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
email: Rainer<at>krugs<dot>de
PGP: 0x0F52F982
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Type: application/pgp-signature
Size: 454 bytes
Desc: not available
URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20151006/0da07742/attachment.bin>


More information about the R-help mailing list