[R] Bug in auglag?

Rainer M Krug Rainer at krugs.de
Tue Oct 6 15:20:42 CEST 2015


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:

--8<---------------cut here---------------start------------->8---
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)
}
--8<---------------cut here---------------end--------------->8---

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, 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:      Rainer at krugs.de

Skype:      RMkrug

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/51e4da9e/attachment.bin>


More information about the R-help mailing list