[R] Bug in auglag?
Ravi Varadhan
ravi.varadhan at jhu.edu
Tue Oct 6 22:58:53 CEST 2015
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 <Rainer at krugs.de>
Sent: Tuesday, October 6, 2015 9:20 AM
To: Ravi Varadhan
Cc: 'r-help at 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:
--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
More information about the R-help
mailing list