The eval() call could also throw an error that would leave the input
environment modified. Better change along the lines described in the
bug report at

https://bugs.r-project.org/bugzilla/show_bug.cgi?id=17831

Best,

luke

On Tue, 16 Jun 2020, Raimundo Neto wrote:

Dear all

As far as I could trace, looking at the function C function numeric_deriv,
this unwanted behavior comes from the inner most loop in, at the very end of
the function,
for(i = 0, start = 0; i < LENGTH(theta); i++) {
  for(j = 0; j < LENGTH(VECTOR_ELT(pars, i)); j++, start += LENGTH(ans)) {
    SEXP ans_del;
    double origPar, xx, delta;

    origPar = REAL(VECTOR_ELT(pars, i))[j];
    xx = fabs(origPar);
    delta = (xx == 0) ? eps : xx*eps;
    REAL(VECTOR_ELT(pars, i))[j] += rDir[i] * delta;
    PROTECT(ans_del = eval(expr, rho));
    if(!isReal(ans_del)) ans_del = coerceVector(ans_del, REALSXP);
    UNPROTECT(1);
    for(k = 0; k < LENGTH(ans); k++) {
      if (!R_FINITE(REAL(ans_del)[k]))
        error(_("Missing value or an infinity produced when evaluating the
model"));
      REAL(gradient)[start + k] = rDir[i] * (REAL(ans_del)[k] -
REAL(ans)[k])/delta;
    }
    REAL(VECTOR_ELT(pars, i))[j] = origPar;
  }
}
Maybe a (naive?) fix is change the if statement in the inner most loop to

if (!R_FINITE(REAL(ans_del)[k])) {
  REAL(VECTOR_ELT(pars, i))[j] = origPar;
  error(_("Missing value or an infinity produced when evaluating the
model"));
}


Regards,
Raimundo Neto


Em ter., 16 de jun. de 2020 às 11:31, <luke-tier...@uiowa.edu> escreveu:
      Thanks; definitely a bug. I've submitted it to the bug tracker
      at

      https://bugs.r-project.org/bugzilla/show_bug.cgi?id=17831

      Best,

      luke

      On Mon, 15 Jun 2020, Raimundo Neto wrote:

      > Dear R developers,
      >
      > I've run into a weird behavior of the numericDeriv function
      (from the stats
      > package) which I also posted on StackOverflow (question has
      same title as
      > this email, except for the version of R).
      >
      > Running the code bellow we can see that the numericDeriv
      function gives an
      > error as the derivative of x^a wrt a is x^a * log(x) and log
      is not defined
      > for negative numbers. However, seems like the function changes
      the value of
      > env1$a from 3 to 3.000000044703483581543. If x is a vector of
      positive
      > values numericDeriv function completes the task without
      errors  and env1$a
      > remains unchanged as expected.
      >
      > This happened to me running R 4.0.1 on Ubuntu 20.04 and also
      to another
      > StackOverflow user using running the same version of R on
      Windows 10. I
      > wonder, is this an intended behavior of the function or really
      a bug?
      >
      > options(digits=22)
      > env1 = new.env()
      > env1$x = rnorm(10)
      > env1$a = 3
      > eval(quote(x^a), env1)
      > numericDeriv(quote(x^a), "a", env1)
      > eval(quote(x^a), env1)
      > env1$a
      >
      > Thank you!
      > Raimundo Neto
      >
      >       [[alternative HTML version deleted]]
      >
      > ______________________________________________
      > R-devel@r-project.org mailing list
      > https://stat.ethz.ch/mailman/listinfo/r-devel
      >

      --
      Luke Tierney
      Ralph E. Wareham Professor of Mathematical Sciences
      University of Iowa                  Phone:           
       319-335-3386
      Department of Statistics and        Fax:             
       319-335-3017
          Actuarial Science
      241 Schaeffer Hall                  email: 
       luke-tier...@uiowa.edu
      Iowa City, IA 52242                 WWW: 
      http://www.stat.uiowa.edu




--
Luke Tierney
Ralph E. Wareham Professor of Mathematical Sciences
University of Iowa                  Phone:             319-335-3386
Department of Statistics and        Fax:               319-335-3017
   Actuarial Science
241 Schaeffer Hall                  email:   luke-tier...@uiowa.edu
Iowa City, IA 52242                 WWW:  http://www.stat.uiowa.edu
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

Reply via email to