Title: | Numerical Tools for 'Rcpp' and Lambda Functions |
---|---|
Description: | Provides a 'C++' API for routinely used numerical tools such as integration, root-finding, and optimization, where function arguments are given as lambdas. This facilitates 'Rcpp' programming, enabling the development of 'R'-like code in 'C++' where functions can be defined on the fly and use variables in the surrounding environment. |
Authors: | Andrew M. Raim [aut, cre] |
Maintainer: | Andrew M. Raim <[email protected]> |
License: | MIT + file LICENSE |
Version: | 0.1.2 |
Built: | 2025-01-26 05:13:23 UTC |
Source: | https://github.com/andrewraim/fntl |
Package documentation
Maintainer: Andrew M. Raim [email protected] (ORCID)
Useful links:
Get an arguments list for internal methods with the default settings. This object can be adjusted and passed to the respective function.
findroot_args() optimize_args() integrate_args() cg_args() bfgs_args() lbfgsb_args() neldermead_args() nlm_args() richardson_args()
findroot_args() optimize_args() integrate_args() cg_args() bfgs_args() lbfgsb_args() neldermead_args() nlm_args() richardson_args()
An argument list corresponding to the specified function. The elements of the list are named and supplied with default values. See the package vignette for further details.
findroot_args
is documented in the section "Root-Finding".
optimize_args
is documented in the section "Univariate Optimization".
integrate_args
is documented in the section "Integration".
cg_args
is documented in the section "Conjugate Gradient".
bfgs_args
is documented in the section "BFGS".
lbfgsb_args
is documented in the section "L-BFGS-B".
neldermead_args
is documented in the section "Nelder-Mead".
nlm_args
is documented in the section "Newton-Type Algorithm for
Nonlinear Optimization".
richardson_args
is documented in the section "Richardson
Extrapolated Finite Differences".
Numerical Derivatives via Finite Differences
fd_deriv1(f, x, i, h, fd_type) fd_deriv2(f, x, i, j, h_i, h_j, fd_type) deriv1(f, x, i, args, fd_type) deriv2(f, x, i, j, args, fd_type)
fd_deriv1(f, x, i, h, fd_type) fd_deriv2(f, x, i, j, h_i, h_j, fd_type) deriv1(f, x, i, args, fd_type) deriv2(f, x, i, j, args, fd_type)
f |
Function to differentiate. |
x |
Scalar at which to evaluate the derivative. |
i |
First coordinate to differentiate. |
h |
Step size in the first coordinate. |
fd_type |
Type of derivative: |
j |
Second coordinate to differentiate. |
h_i |
Step size in the first coordinate. |
h_j |
Step size in the second coordinate. |
args |
List of additional arguments from the function |
fd_deriv1
and fd_deriv2
return a single numeric value corresponding to
the first and second derivative via finite differences. deriv1
and
deriv2
return a list with the form of a richardson_result
described in
section "Richardson Extrapolated Finite Differences" of the package
vignette.
args = richardson_args() f = sin # Try 2nd derivatives of a univariate function x0 = 0.5 print(-sin(x0)) ## Exact answer for f''(x0) fd_deriv2(f, x0, i = 0, j = 0, h_i = 0.001, h_j = 0.001, fd_type = 0) fd_deriv2(f, x0, i = 0, j = 0, h_i = 0.001, h_j = 0.001, fd_type = 1) fd_deriv2(f, x0, i = 0, j = 0, h_i = 0.001, h_j = 0.001, fd_type = 2) deriv2(f, x0, i = 0, j = 0, args, fd_type = 0) # Try 2nd derivatives of a bivariate function f = function(x) { sin(x[1]) + cos(x[2]) } x0 = c(0.5, 0.25) print(-sin(x0[1])) ## Exact answer for f_xx(x0) print(-cos(x0[2])) ## Exact answer for f_yy(x0) print(0) ## Exact answer for f_xy(x0,y0) numDeriv::hessian(f, x0) fd_deriv2(f, x0, i = 0, j = 0, h_i = 0.001, h_j = 0.001, fd_type = 0) fd_deriv2(f, x0, i = 0, j = 0, h_i = 0.001, h_j = 0.001, fd_type = 1) fd_deriv2(f, x0, i = 0, j = 0, h_i = 0.001, h_j = 0.001, fd_type = 2) fd_deriv2(f, x0, i = 0, j = 1, h_i = 0.001, h_j = 0.001, fd_type = 0) fd_deriv2(f, x0, i = 0, j = 1, h_i = 0.001, h_j = 0.001, fd_type = 1) fd_deriv2(f, x0, i = 0, j = 1, h_i = 0.001, h_j = 0.001, fd_type = 2) fd_deriv2(f, x0, i = 1, j = 1, h_i = 0.001, h_j = 0.001, fd_type = 0) fd_deriv2(f, x0, i = 1, j = 1, h_i = 0.001, h_j = 0.001, fd_type = 1) fd_deriv2(f, x0, i = 1, j = 1, h_i = 0.001, h_j = 0.001, fd_type = 2) deriv2(f, x0, i = 1, j = 1, args, fd_type = 0) deriv2(f, x0, i = 1, j = 1, args, fd_type = 1) deriv2(f, x0, i = 1, j = 1, args, fd_type = 2)
args = richardson_args() f = sin # Try 2nd derivatives of a univariate function x0 = 0.5 print(-sin(x0)) ## Exact answer for f''(x0) fd_deriv2(f, x0, i = 0, j = 0, h_i = 0.001, h_j = 0.001, fd_type = 0) fd_deriv2(f, x0, i = 0, j = 0, h_i = 0.001, h_j = 0.001, fd_type = 1) fd_deriv2(f, x0, i = 0, j = 0, h_i = 0.001, h_j = 0.001, fd_type = 2) deriv2(f, x0, i = 0, j = 0, args, fd_type = 0) # Try 2nd derivatives of a bivariate function f = function(x) { sin(x[1]) + cos(x[2]) } x0 = c(0.5, 0.25) print(-sin(x0[1])) ## Exact answer for f_xx(x0) print(-cos(x0[2])) ## Exact answer for f_yy(x0) print(0) ## Exact answer for f_xy(x0,y0) numDeriv::hessian(f, x0) fd_deriv2(f, x0, i = 0, j = 0, h_i = 0.001, h_j = 0.001, fd_type = 0) fd_deriv2(f, x0, i = 0, j = 0, h_i = 0.001, h_j = 0.001, fd_type = 1) fd_deriv2(f, x0, i = 0, j = 0, h_i = 0.001, h_j = 0.001, fd_type = 2) fd_deriv2(f, x0, i = 0, j = 1, h_i = 0.001, h_j = 0.001, fd_type = 0) fd_deriv2(f, x0, i = 0, j = 1, h_i = 0.001, h_j = 0.001, fd_type = 1) fd_deriv2(f, x0, i = 0, j = 1, h_i = 0.001, h_j = 0.001, fd_type = 2) fd_deriv2(f, x0, i = 1, j = 1, h_i = 0.001, h_j = 0.001, fd_type = 0) fd_deriv2(f, x0, i = 1, j = 1, h_i = 0.001, h_j = 0.001, fd_type = 1) fd_deriv2(f, x0, i = 1, j = 1, h_i = 0.001, h_j = 0.001, fd_type = 2) deriv2(f, x0, i = 1, j = 1, args, fd_type = 0) deriv2(f, x0, i = 1, j = 1, args, fd_type = 1) deriv2(f, x0, i = 1, j = 1, args, fd_type = 2)
Find Root
findroot_bisect(f, lower, upper, args) findroot_brent(f, lower, upper, args)
findroot_bisect(f, lower, upper, args) findroot_brent(f, lower, upper, args)
f |
Function for which a root is desired. |
lower |
Lower limit of search interval. Must be finite. |
upper |
Upper limit of search interval. Must be finite. |
args |
List of additional arguments from the function |
A list with the form of a findroot_result
described in section
"Root-Finding" of the package vignette.
f = function(x) { x^2 - 1 } args = findroot_args() findroot_bisect(f, 0, 10, args) findroot_brent(f, 0, 10, args)
f = function(x) { x^2 - 1 } args = findroot_args() findroot_bisect(f, 0, 10, args) findroot_brent(f, 0, 10, args)
Numerical Gradient Vector
gradient0(f, x, args)
gradient0(f, x, args)
f |
Function to differentiate. |
x |
Vector at which to evaluate the gradient. |
args |
List of additional arguments from the function |
A list with the form of a gradient_result
described in section "Gradient"
of the package vignette.
f = function(x) { sum(sin(x)) } args = richardson_args() x0 = seq(0, 1, length.out = 5) cos(x0) ## Exact answer gradient0(f, x0, args) numDeriv::grad(f, x0)
f = function(x) { sum(sin(x)) } args = richardson_args() x0 = seq(0, 1, length.out = 5) cos(x0) ## Exact answer gradient0(f, x0, args) numDeriv::grad(f, x0)
Numerical Hessian
hessian0(f, x, args)
hessian0(f, x, args)
f |
Function to differentiate. |
x |
Vector at which to evaluate the Hessian. |
args |
List of additional arguments from the function |
A list with the form of a hessian_result
described in section "Hessian"
of the package vignette.
f = function(x) { sum(x^2) } x0 = seq(1, 10, length.out = 5) args = richardson_args() hessian0(f, x0, args) numDeriv::hessian(f, x0)
f = function(x) { sum(x^2) } x0 = seq(1, 10, length.out = 5) args = richardson_args() hessian0(f, x0, args) numDeriv::hessian(f, x0)
Compute the integral .
integrate0(f, lower, upper, args)
integrate0(f, lower, upper, args)
f |
Function to integrate. |
lower |
Lower limit of integral. |
upper |
Upper limit of integral. |
args |
List of additional arguments from the function |
A list with the form of a integrate_result
described in section
"Integration" of the package vignette.
f = function(x) { exp(-x^2 / 2) } args = integrate_args() integrate0(f, 0, 10, args)
f = function(x) { exp(-x^2 / 2) } args = integrate_args() integrate0(f, 0, 10, args)
Numerical Jacobian Matrix
jacobian0(f, x, args)
jacobian0(f, x, args)
f |
Function to differentiate. |
x |
Vector at which to evaluate the Jacobian. |
args |
List of additional arguments from the function |
A list with the form of a jacobian_result
described in section "Jacobian"
of the package vignette.
f = function(x) { cumsum(sin(x)) } x0 = seq(1, 10, length.out = 5) args = richardson_args() out = jacobian0(f, x0, args) print(out$value) numDeriv::jacobian(f, x0)
f = function(x) { cumsum(sin(x)) } x0 = seq(1, 10, length.out = 5) args = richardson_args() out = jacobian0(f, x0, args) print(out$value) numDeriv::jacobian(f, x0)
Matrix Apply Functions
mat_apply(X, f) row_apply(X, f) col_apply(X, f)
mat_apply(X, f) row_apply(X, f) col_apply(X, f)
X |
A matrix |
f |
The function to apply. |
The mat_apply
, row_apply
, and col_apply
C++ functions are intended to
operate like the following calls in R, respectively.
apply(x, c(1,2), f) apply(x, 1, f) apply(x, 2, f)
The R functions exposed here are specific to numeric-valued matrices, but the underlying C++ functions are intended to work with any type of Rcpp Matrix.
mat_apply
returns a matrix. row_apply
and col_apply
return a vector.
See section "Apply" of the package vignette for details.
X = matrix(1:12, nrow = 4, ncol = 3) mat_apply(X, f = function(x) { x^(1/3) }) row_apply(X, f = function(x) { sum(x^2) }) col_apply(X, f = function(x) { sum(x^2) })
X = matrix(1:12, nrow = 4, ncol = 3) mat_apply(X, f = function(x) { x^(1/3) }) row_apply(X, f = function(x) { sum(x^2) }) col_apply(X, f = function(x) { sum(x^2) })
Multivariate Optimization
cg1(init, f, g, args) cg2(init, f, args) bfgs1(init, f, g, args) bfgs2(init, f, args) lbfgsb1(init, f, g, args) lbfgsb2(init, f, args) neldermead(init, f, args) nlm1(init, f, g, h, args) nlm2(init, f, g, args) nlm3(init, f, args)
cg1(init, f, g, args) cg2(init, f, args) bfgs1(init, f, g, args) bfgs2(init, f, args) lbfgsb1(init, f, g, args) lbfgsb2(init, f, args) neldermead(init, f, args) nlm1(init, f, g, h, args) nlm2(init, f, g, args) nlm3(init, f, args)
init |
Initial value |
f |
Function |
g |
Gradient function of |
args |
List of additional arguments for optimization. |
h |
Hessian function of |
The argument args
should be a list constructed from one of the following
functions:
bfgs_args
for BFGS;
lbfgsb_args
for L-BFGS-B;
cg_args
for CG;
neldermead_args
for Nelder-Mead;
nlm_args
for the Newton-type algorithm used in nlm
.
When g
or h
are omitted, the gradient or Hessian will be respectively
be computed via finite differences.
A list with results corresponding to the specified function. See the package vignette for further details.
cg1
and cg2
return a cg_result
which is documented in the section
"Conjugate Gradient".
bfgs1
and bfgs2
return a bfgs_result
which is documented in the
section "BFGS".
lbfgsb1
and lbfgsb2
return a lbfgsb_result
which is documented in
the section "L-BFGS-B".
neldermead
returns a neldermead_result
which is documented in
the section "Nelder-Mead".
nlm1
, nlm2
, and nlm3
return a nlm_result
which is documented in
the section "Newton-Type Algorithm for Nonlinear Optimization".
f = function(x) { sum(x^2) } g = function(x) { 2*x } h = function(x) { 2*diag(length(x)) } x0 = c(1,1) args = cg_args() cg1(x0, f, g, args) cg2(x0, f, args) args = bfgs_args() bfgs1(x0, f, g, args) bfgs2(x0, f, args) args = lbfgsb_args() lbfgsb1(x0, f, g, args) lbfgsb2(x0, f, args) args = neldermead_args() neldermead(x0, f, args) args = nlm_args() nlm1(x0, f, g, h, args) nlm2(x0, f, g, args) nlm3(x0, f, args)
f = function(x) { sum(x^2) } g = function(x) { 2*x } h = function(x) { 2*diag(length(x)) } x0 = c(1,1) args = cg_args() cg1(x0, f, g, args) cg2(x0, f, args) args = bfgs_args() bfgs1(x0, f, g, args) bfgs2(x0, f, args) args = lbfgsb_args() lbfgsb1(x0, f, g, args) lbfgsb2(x0, f, args) args = neldermead_args() neldermead(x0, f, args) args = nlm_args() nlm1(x0, f, g, h, args) nlm2(x0, f, g, args) nlm3(x0, f, args)
Compute "outer" matrices and matrix-vector products based on a function that operators on pairs of rows. See details.
outer1(X, f) outer2(X, Y, f) outer1_matvec(X, f, a) outer2_matvec(X, Y, f, a)
outer1(X, f) outer2(X, Y, f) outer1_matvec(X, f, a) outer2_matvec(X, Y, f, a)
X |
A numerical matrix. |
f |
Function |
Y |
A numerical matrix. |
a |
A scalar vector. |
The outer1
function computes the symmetric matrix
and the outer1_matvec
operation computes the -dimensional vector
The outer2
operation computes the matrix
and the outer2_matvec
operation computes the -dimensional vector
outer1
and outer2
return a matrix. outer1_matvec
and outer2_matvec
return a vector. See section "Outer" of the package vignette for details.
set.seed(1234) f = function(x,y) { sum( (x - y)^2 ) } X = matrix(rnorm(12), 6, 2) Y = matrix(rnorm(10), 5, 2) outer1(X, f) outer2(X, Y, f) a = rep(1, 6) b = rep(1, 5) outer1_matvec(X, f, a) outer2_matvec(X, Y, f, b)
set.seed(1234) f = function(x,y) { sum( (x - y)^2 ) } X = matrix(rnorm(12), 6, 2) Y = matrix(rnorm(10), 5, 2) outer1(X, f) outer2(X, Y, f) a = rep(1, 6) b = rep(1, 5) outer1_matvec(X, f, a) outer2_matvec(X, Y, f, b)
Solve the system where
is a matrix-free
representation of the linear operation
.
solve_cg(l, b, init, args)
solve_cg(l, b, init, args)
l |
A linear transformation of |
b |
A vector. |
init |
Initial value of solution. |
args |
List of additional arguments from |
A list with the form of a solve_cg_result
described in section "Conjugate
Gradient" of the package vignette.
set.seed(1234) n = 8 idx_diag = cbind(1:n, 1:n) idx_ldiag = cbind(2:n, 1:(n-1)) idx_udiag = cbind(1:(n-1), 2:n) b = rep(1, n) ## Solution by explicit computation of solve(A, b) A = matrix(0, n, n) A[idx_diag] = 2 A[idx_ldiag] = 1 A[idx_udiag] = 1 solve(A, b) ## Solve iteratively with solve_cg f = function(x) { A %*% x } args = cg_args() init = rep(0, n) solve_cg(f, b, init, args)
set.seed(1234) n = 8 idx_diag = cbind(1:n, 1:n) idx_ldiag = cbind(2:n, 1:(n-1)) idx_udiag = cbind(1:(n-1), 2:n) b = rep(1, n) ## Solution by explicit computation of solve(A, b) A = matrix(0, n, n) A[idx_diag] = 2 A[idx_ldiag] = 1 A[idx_udiag] = 1 solve(A, b) ## Solve iteratively with solve_cg f = function(x) { A %*% x } args = cg_args() init = rep(0, n) solve_cg(f, b, init, args)
Density, CDF, quantile, and drawing functions for a univariate distribution
with density , cdf,
, and quantile function
truncated to the interval
.
d_trunc(x, lo, hi, f, F, log = FALSE) p_trunc(x, lo, hi, F, lower = TRUE, log = FALSE) q_trunc(p, lo, hi, F, Finv, lower = TRUE, log = FALSE) r_trunc(n, lo, hi, F, Finv)
d_trunc(x, lo, hi, f, F, log = FALSE) p_trunc(x, lo, hi, F, lower = TRUE, log = FALSE) q_trunc(p, lo, hi, F, Finv, lower = TRUE, log = FALSE) r_trunc(n, lo, hi, F, Finv)
x |
Vector of quantiles. |
lo |
Vector of lower limits. |
hi |
Vector of upper limits. |
f |
Density function with form |
F |
CDF function with signature |
log |
logical; if |
lower |
logical; if |
p |
Vector of probabilities. |
Finv |
Quantile function with signature |
n |
Number of draws. |
Vector with results.
d_trunc
computes the density, r_trunc
generates random deviates,
p_trunc
computes the CDF, and q_trunc
computes quantiles.
library(tidyverse) m = 100 ## Length of sequence for density, CDF, etc shape1 = 5 shape2 = 2 lo = 0.5 hi = 0.7 # Density, CDF, and quantile function for untruncated distribution ff = function(x, log) { dbeta(x, shape1, shape2, log = log) } FF = function(x, lower, log) { pbeta(x, shape1, shape2, lower.tail = lower, log = log) } FFinv = function(x, lower, log) { qbeta(x, shape1, shape2, lower.tail = lower, log = log) } # Compare truncated and untruncated densities xseq = seq(0, 1, length.out = m) lo_vec = rep(lo, m) hi_vec = rep(hi, m) f0seq = ff(xseq, log = FALSE) fseq = d_trunc(xseq, lo_vec, hi_vec, ff, FF) data.frame(x = xseq, f = fseq, f0 = f0seq) %>% ggplot() + geom_line(aes(xseq, fseq)) + geom_line(aes(xseq, f0seq), lty = 2) + xlab("x") + ylab("Density") + theme_minimal() # Compare truncated densities and empirical density of generated draws n = 100000 lo_vec = rep(lo, n) hi_vec = rep(hi, n) x = r_trunc(n = n, lo_vec, hi_vec, FF, FFinv) hist(x, probability = TRUE, breaks = 15) points(xseq, fseq) # Compare empirical CDF of draws with CDF function Femp = ecdf(x) lo_vec = rep(lo, m) hi_vec = rep(hi, m) Fseq = p_trunc(xseq, lo_vec, hi_vec, FF) data.frame(x = xseq, FF = Fseq) %>% mutate(F0 = Femp(x)) %>% ggplot() + geom_line(aes(xseq, FF), lwd = 1.2) + geom_line(aes(xseq, F0), col = "orange") + xlab("x") + ylab("Probability") + theme_minimal() # Compare empirical quantiles of draws with quantile function pseq = seq(0, 1, length.out = m) lo_vec = rep(lo, m) hi_vec = rep(hi, m) Finvseq = q_trunc(pseq, lo_vec, hi_vec, FF, FFinv) Finvemp = quantile(x, prob = pseq) data.frame(p = pseq, Finv = Finvseq, Finvemp = Finvemp) %>% ggplot() + geom_line(aes(pseq, Finv), lwd = 1.2) + geom_line(aes(pseq, Finvemp), col = "orange") + xlab("p") + ylab("Quantile") + theme_minimal()
library(tidyverse) m = 100 ## Length of sequence for density, CDF, etc shape1 = 5 shape2 = 2 lo = 0.5 hi = 0.7 # Density, CDF, and quantile function for untruncated distribution ff = function(x, log) { dbeta(x, shape1, shape2, log = log) } FF = function(x, lower, log) { pbeta(x, shape1, shape2, lower.tail = lower, log = log) } FFinv = function(x, lower, log) { qbeta(x, shape1, shape2, lower.tail = lower, log = log) } # Compare truncated and untruncated densities xseq = seq(0, 1, length.out = m) lo_vec = rep(lo, m) hi_vec = rep(hi, m) f0seq = ff(xseq, log = FALSE) fseq = d_trunc(xseq, lo_vec, hi_vec, ff, FF) data.frame(x = xseq, f = fseq, f0 = f0seq) %>% ggplot() + geom_line(aes(xseq, fseq)) + geom_line(aes(xseq, f0seq), lty = 2) + xlab("x") + ylab("Density") + theme_minimal() # Compare truncated densities and empirical density of generated draws n = 100000 lo_vec = rep(lo, n) hi_vec = rep(hi, n) x = r_trunc(n = n, lo_vec, hi_vec, FF, FFinv) hist(x, probability = TRUE, breaks = 15) points(xseq, fseq) # Compare empirical CDF of draws with CDF function Femp = ecdf(x) lo_vec = rep(lo, m) hi_vec = rep(hi, m) Fseq = p_trunc(xseq, lo_vec, hi_vec, FF) data.frame(x = xseq, FF = Fseq) %>% mutate(F0 = Femp(x)) %>% ggplot() + geom_line(aes(xseq, FF), lwd = 1.2) + geom_line(aes(xseq, F0), col = "orange") + xlab("x") + ylab("Probability") + theme_minimal() # Compare empirical quantiles of draws with quantile function pseq = seq(0, 1, length.out = m) lo_vec = rep(lo, m) hi_vec = rep(hi, m) Finvseq = q_trunc(pseq, lo_vec, hi_vec, FF, FFinv) Finvemp = quantile(x, prob = pseq) data.frame(p = pseq, Finv = Finvseq, Finvemp = Finvemp) %>% ggplot() + geom_line(aes(pseq, Finv), lwd = 1.2) + geom_line(aes(pseq, Finvemp), col = "orange") + xlab("p") + ylab("Quantile") + theme_minimal()
Univariate Optimization
goldensection(f, lower, upper, args) optimize_brent(f, lower, upper, args)
goldensection(f, lower, upper, args) optimize_brent(f, lower, upper, args)
f |
Function to optimize. |
lower |
Lower limit of search interval. Must be finite. |
upper |
Upper limit of search interval. Must be finite. |
args |
List of additional arguments from the function |
A list with the form of a optimize_result
described in section
"Univariate Optimization" of the package vignette.
f = function(x) { x^2 - 1 } args = optimize_args() goldensection(f, 0, 10, args) optimize_brent(f, 0, 10, args)
f = function(x) { x^2 - 1 } args = optimize_args() goldensection(f, 0, 10, args) optimize_brent(f, 0, 10, args)
Matrix Which Function
which0(X, f)
which0(X, f)
X |
A matrix |
f |
A predicate to apply to each element of |
The which
C++ functions are intended to operate like the following call
in R.
which(f(X), arr.ind = TRUE) - 1
The R functions exposed here are specific to numeric-valued matrices, but the underlying C++ functions are intended to work with any type of Rcpp Matrix.
A matrix with two columns. Each row contains a row and column index
corresponding to an element of that matches the criteria of
.
See section "Which" of the package vignette for details.
X = matrix(1:12 / 6, nrow = 4, ncol = 3) f = function(x) { x < 1 } which0(X, f)
X = matrix(1:12 / 6, nrow = 4, ncol = 3) f = function(x) { x < 1 } which0(X, f)