Commit 561acd6c authored by Radford Neal's avatar Radford Neal

add deriv extensions from R-3.4.0, except for sinpi, etc. Also stop

putting in parens in deriv.
parent 62bce3e8
......@@ -6,6 +6,25 @@
\encoding{UTF-8}
\section{CHANGES IN VERSION RELEASED 2018-00-00}{
\subsection{FEATURE CHANGES}{
\itemize{
\item \code{D} and \code{deriv} no longer add parenthesis in (some of)
the places where they would be needed according to precedence, since
this slows evaluation, and is unnecessary considering that
\code{deparse} adds such parentheses.
}}
\subsection{NEW FEATURES FROM R CORE RELEASES}{
\itemize{
\item From R-3.4.0: The \code{deriv()} and similar functions now can compute
derivatives of \code{log1p()}, \code{sinpi()} and similar
one-argument functions, thanks to a contribution by Jerry Lewis.
[ Except that \code{sinpi}, etc. are not yet in pqR. ]
}}
}
\section{CHANGES IN VERSION RELEASED 2018-11-18}{
\subsection{INTRODUCTION}{
......
Incorporate extensions to deriv from R-3.4.0.
% File src/library/stats/man/deriv.Rd
% Part of the R package, http://www.R-project.org
% Copyright 1995-2007 R Core Team
% Modifications for pqR Copyright (c) 2018 Radford M. Neal.
% Distributed under GPL 2 or later
\name{deriv}
......@@ -66,9 +67,11 @@ deriv3(expr, \dots)
The internal code knows about the arithmetic operators \code{+},
\code{-}, \code{*}, \code{/} and \code{^}, and the single-variable
functions \code{exp}, \code{log}, \code{sin}, \code{cos}, \code{tan},
functions \code{exp}, \code{expm1}, \code{log}, \code{log2}, \code{log10},
\code{log1p}, \code{sin}, \code{cos}, \code{tan},
\code{sinh}, \code{cosh}, \code{sqrt}, \code{pnorm}, \code{dnorm},
\code{asin}, \code{acos}, \code{atan}, \code{gamma}, \code{lgamma},
\code{asin}, \code{acos}, \code{atan}, \code{factorial}, \code{lfactorial},
\code{gamma}, \code{lgamma},
\code{digamma} and \code{trigamma}, as well as \code{psigamma} for one
or two arguments (but derivative only with respect to the first).
(Note that only the standard normal distribution is considered.)
......@@ -145,11 +148,25 @@ DD <- function(expr,name, order = 1) {
else DD(D(expr, name), name, order - 1)
}
DD(expression(sin(x^2)), "x", 3)
## showing the limits of the internal "simplify()" :
\dontrun{
-sin(x^2) * (2 * x) * 2 + ((cos(x^2) * (2 * x) * (2 * x) + sin(x^2) *
\dontrun{-sin(x^2) * (2 * x) * 2 + ((cos(x^2) * (2 * x) * (2 * x) + sin(x^2) *
2) * (2 * x) + sin(x^2) * (2 * x) * 2)
}
D(quote(log1p(x^2)), "x") ## log1p(x) = log(1 + x)
stopifnot(identical(
D(quote(log1p(x^2)), "x"),
D(quote(log(1+x^2)), "x")))
D(quote(expm1(x^2)), "x") ## expm1(x) = exp(x) - 1
stopifnot(identical(
D(quote(expm1(x^2)), "x") -> Dex1,
D(quote(exp(x^2)-1), "x")),
identical(Dex1, quote(exp(x^2) * (2 * x))))
stopifnot(identical(D(quote(log2 (x^2)), "x"),
quote(2 * x/(x^2 * log(2)))),
identical(D(quote(log10(x^2)), "x"),
quote(2 * x/(x^2 * log(10)))))
}
\keyword{math}
\keyword{nonlinear}
......
/*
* pqR : A pretty quick version of R
* Copyright (C) 2013, 2014 by Radford M. Neal
* Copyright (C) 2013, 2014, 2018 by Radford M. Neal
*
* Based on R : A Computer Language for Statistical Data Analysis
* Copyright (C) 1995, 1996 Robert Gentleman and Ross Ihaka
* Copyright (C) 1998-2011 The R Core Team.
* Copyright (C) 2004-5 The R Foundation
*
* Some extensions from R-3.4.0 are incorporated,
* Copyright (C) 2011-2017 The R Core Team
*
* The changes in pqR from R-2.15.0 distributed by the R Core Team are
* documented in the NEWS and MODS files in the top-level source directory.
*
......@@ -61,6 +64,14 @@ static SEXP LGammaSymbol;
static SEXP DiGammaSymbol;
static SEXP TriGammaSymbol;
static SEXP PsiSymbol;
static SEXP ExpM1Symbol;
static SEXP Log1PSymbol;
static SEXP Log2Symbol;
static SEXP Log10Symbol;
static SEXP FactorialSymbol;
static SEXP LFactorialSymbol;
/* sinpi, etc. from R-3.4.0 not yet included,
since these functions aren't yet in pqR. */
static Rboolean Initialized = FALSE;
......@@ -94,6 +105,12 @@ static void InitDerivSymbols(void)
DiGammaSymbol = install("digamma");
TriGammaSymbol = install("trigamma");
PsiSymbol = install("psigamma");
ExpM1Symbol = install("expm1");
Log1PSymbol = install("log1p");
Log2Symbol = install("log2");
Log10Symbol = install("log10");
FactorialSymbol = install("factorial");
LFactorialSymbol = install("lfactorial");
Initialized = TRUE;
}
......@@ -262,6 +279,18 @@ static SEXP simplify(SEXP fun, SEXP arg1, SEXP arg2)
if (arg2 == R_MissingArg) ans = lang2(PsiSymbol, arg1);
else ans = lang3(PsiSymbol, arg1, arg2);
}
else if (fun == ExpM1Symbol) {
/* FIXME: simplify expm1(log1p( E )) = E */
ans = lang2(ExpM1Symbol, arg1);
}
else if (fun == Log1PSymbol) {
/* FIXME: simplify log1p(expm1( E )) = E */
ans = lang2(Log1PSymbol, arg1);
}
else if (fun == Log2Symbol) ans = lang2(Log2Symbol, arg1);
else if (fun == Log10Symbol) ans = lang2(Log10Symbol, arg1);
else if (fun == FactorialSymbol)ans = lang2(FactorialSymbol, arg1);
else if (fun == LFactorialSymbol)ans = lang2(LFactorialSymbol, arg1);
else ans = Constant(NA_REAL);
/* FIXME */
#ifdef NOTYET
......@@ -527,7 +556,50 @@ static SEXP D(SEXP expr, SEXP var)
UNPROTECT(3);
}
}
else if (CAR(expr) == ExpM1Symbol) {
ans = simplify(TimesSymbol,
PP_S2(ExpSymbol, CADR(expr)),
PP(D(CADR(expr), var)));
UNPROTECT(2);
}
else if (CAR(expr) == Log1PSymbol) {
ans = simplify(DivideSymbol,
PP(D(CADR(expr), var)),
PP_S(PlusSymbol, PP(Constant(1.)), CADR(expr)));
UNPROTECT(3);
}
else if (CAR(expr) == Log2Symbol) {
ans = simplify(DivideSymbol,
PP(D(CADR(expr), var)),
PP_S(TimesSymbol, CADR(expr),
PP_S2(LogSymbol, PP(Constant(2.)))));
UNPROTECT(4);
}
else if (CAR(expr) == Log10Symbol) {
ans = simplify(DivideSymbol,
PP(D(CADR(expr), var)),
PP_S(TimesSymbol, CADR(expr),
PP_S2(LogSymbol, PP(Constant(10.)))));
UNPROTECT(4);
}
else if (CAR(expr) == LFactorialSymbol) {
ans = simplify(TimesSymbol,
PP(D(CADR(expr), var)),
PP_S2(DiGammaSymbol, PP_S(PlusSymbol,
CADR(expr),
PP(ScalarInteger(1)))));
UNPROTECT(4);
}
else if (CAR(expr) == FactorialSymbol) {
ans = simplify(TimesSymbol,
PP(D(CADR(expr), var)),
PP_S(TimesSymbol,
expr,
PP_S2(DiGammaSymbol, PP_S(PlusSymbol,
CADR(expr),
PP(ScalarInteger(1))))));
UNPROTECT(5);
}
else {
SEXP u = deparse1(CAR(expr), 0, SIMPLEDEPARSE);
error(_("Function '%s' is not in the derivatives table"),
......@@ -579,14 +651,15 @@ static int isPowerForm(SEXP expr)
&& CAR(expr) == PowerSymbol);
}
/* Add parentheses is some places where they are needed to get the right
precedence. But not all. It's not clear this is needed, since
deparsing puts parentheses in where needed, and evaluation gets it
right with or without parentheses. But we'll do it anyway, as before,
but without altering the expression passed, as it used to do... */
/* Used to add parentheses is some places where they are needed to get
the right precedence. But not all. Shouldn't be needed, since
deparsing puts parentheses in where needed, and evaluation gets it
right with or without parentheses. So this is now disabled. */
static SEXP AddParens(SEXP expr)
{
#if 0 /* NOW DISABLED */
if (TYPEOF(expr) != LANGSXP)
return expr;
......@@ -638,6 +711,9 @@ static SEXP AddParens(SEXP expr)
}
}
UNPROTECT(1);
#endif
return expr;
}
......
......@@ -23,7 +23,7 @@ Type 'demo()' for some demos, 'help()' for on-line help, or
Type 'q()' to quit R.
3 helper threads, task merging enabled, uncompressed pointers.
No helper threads, task merging enabled, uncompressed pointers.
> pkgname <- "stats"
> source(file.path(R.home("share"), "R", "examples-header.R"))
......@@ -3708,7 +3708,7 @@ detaching 'package:cluster'
>
> asOneSidedFormula("age")
~age
<environment: 55c963b15400>
<environment: 56256b2681c0>
> asOneSidedFormula(~ age)
~age
>
......@@ -5978,11 +5978,28 @@ function (b0, b1, th, x)
> DD(expression(sin(x^2)), "x", 3)
-(sin(x^2) * (2 * x) * 2 + ((cos(x^2) * (2 * x) * (2 * x) + sin(x^2) *
2) * (2 * x) + sin(x^2) * (2 * x) * 2))
>
> ## showing the limits of the internal "simplify()" :
> ## Not run:
> ##D -sin(x^2) * (2 * x) * 2 + ((cos(x^2) * (2 * x) * (2 * x) + sin(x^2) *
> ##D 2) * (2 * x) + sin(x^2) * (2 * x) * 2)
> ## End(Not run)
> D(quote(log1p(x^2)), "x") ## log1p(x) = log(1 + x)
2 * x/(1 + x^2)
> stopifnot(identical(
+ D(quote(log1p(x^2)), "x"),
+ D(quote(log(1+x^2)), "x")))
> D(quote(expm1(x^2)), "x") ## expm1(x) = exp(x) - 1
exp(x^2) * (2 * x)
> stopifnot(identical(
+ D(quote(expm1(x^2)), "x") -> Dex1,
+ D(quote(exp(x^2)-1), "x")),
+ identical(Dex1, quote(exp(x^2) * (2 * x))))
>
> stopifnot(identical(D(quote(log2 (x^2)), "x"),
+ quote(2 * x/(x^2 * log(2)))),
+ identical(D(quote(log10(x^2)), "x"),
+ quote(2 * x/(x^2 * log(10)))))
>
>
>
......@@ -6695,7 +6712,7 @@ List of 12
> gf$linkinv
function (eta)
1/eta
<environment: 55c96905cac0>
<environment: 56257091ba80>
> gf$variance(-3:4) #- == (.)^2
[1] 9 4 1 0 1 4 9 16
>
......@@ -7097,7 +7114,7 @@ attr(,".Environment")
> environment(as.formula("y ~ x"))
<environment: R_GlobalEnv>
> environment(as.formula("y ~ x", env=new.env()))
<environment: 55c96a425240>
<environment: 562571b7e840>
>
>
> ## Create a formula for a model with a large number of variables:
......@@ -11511,22 +11528,22 @@ attr(,"degree")
$linkfun
function (mu)
mu
<environment: 55c963d3dc40>
<environment: 56256b659fc0>
$linkinv
function (eta)
eta
<environment: 55c963d3dc40>
<environment: 56256b659fc0>
$mu.eta
function (eta)
rep.int(1, length(eta))
<environment: 55c963d3dc40>
<environment: 56256b659fc0>
$valideta
function (eta)
TRUE
<environment: 55c963d3dc40>
<environment: 56256b659fc0>
$name
[1] "identity"
......@@ -11537,12 +11554,12 @@ attr(,"class")
$linkfun
function (mu)
mu^lambda
<environment: 55c963d3b0c0>
<environment: 56256b657040>
$linkinv
function (eta)
pmax(eta^(1/lambda), .Machine$double.eps)
<environment: 55c963d3b0c0>
<environment: 56256b657040>
>
>
......@@ -15207,7 +15224,7 @@ function (v)
xout = as.double(v), as.integer(length(v)), as.integer(method),
as.double(yleft), as.double(yright), as.double(f), NAOK = TRUE,
PACKAGE = "stats")$xout
<environment: 55c96956a040>
<environment: 562570d7cb40>
attr(,"call")
stepfun(1:3, y0, f = 0)
> ls(envir = environment(sfun0))
......@@ -17423,7 +17440,7 @@ Number of Fisher Scoring iterations: 6
> ### * <FOOTER>
> ###
> cat("Time elapsed: ", proc.time() - get("ptime", pos = 'CheckExEnv'),"\n")
Time elapsed: 6.287 0.111 2.535 0 0
Time elapsed: 2.403 0.1 2.504 0 0
> grDevices::dev.off()
null device
1
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment