Suppose I want to construct the following function:
f <- function(beta) c(y[1]*beta[1]+z[1]*1/beta[2],
y[2]*beta[1]+z[2]*1/beta[2],
: : : :
y[i]*beta[1]^2+z[i]*1/beta[2])
Suppose I have the following data.
y = 1:10
z = 10:19
f <- function(beta) cbind(y) %*% beta^2
jacobian(f, c(1)) #where c(1) is the value for beta.
g <- function(beta) cbind(z) %*% 1/beta
jacobian(g, c(1)) #where c(1) is the value for beta.
Which yields the desired outputs for f and g respectively:
[,1]
[1,] 2
[2,] 4
[3,] 6
[4,] 8
[5,] 10
[6,] 12
[7,] 14
[8,] 16
[9,] 18
[10,] 20
#and
[,1]
[1,] -10
[2,] -11
[3,] -12
[4,] -13
[5,] -14
[6,] -15
[7,] -16
[8,] -17
[9,] -18
[10,] -19
Now I could just merge these two matrices to obtain the jacobian of f and g. However, I would like just one function for the desired output.
I have tried the following but this does not yield the outcome I would like:
u <- function(beta) (cbind(y, z) %*% cbind(beta^2,1/beta))
jacobian(u, c(1,1))
Gives the incorrect output:
[,1] [,2]
[1,] 2 20
[2,] 4 22
[3,] 6 24
[4,] 8 26
[5,] 10 28
[6,] 12 30
[7,] 14 32
[8,] 16 34
[9,] 18 36
[10,] 20 38
[11,] -1 -10
[12,] -2 -11
[13,] -3 -12
[14,] -4 -13
[15,] -5 -14
[16,] -6 -15
[17,] -7 -16
[18,] -8 -17
[19,] -9 -18
[20,] -10 -19
Does anyone know how I can combine function f and g such that I get a 10 x 2 Jacobian matrix?
The jacobian function is structured as follow
library('pracma')
jacobian(f, x0, heps = .Machine$double.eps^(1/3), ...)
f: m functions of n variables.
x0: Numeric vector of length n.
heps: This is h in the derivative formula.
jacobian(): Computes the derivative of each function f_j by variable x_i separately, taking the discrete step h.
The desired output I would like to obtain is
[,1] [,2]
[1,] 2 -10
[2,] 4 -11
[3,] 6 -12
[4,] 8 -13
[5,] 10 -14
[6,] 12 -15
[7,] 14 -16
[8,] 16 -17
[9,] 18 -18
[10,] 20 -19
You went wrong in one particular place:
u <- function(beta) (cbind(y, z) %*% cbind(beta^2,1/beta)) # ^^^^^^^^^^^^^^^^^^^^ # HERE
You used cbind(beta^2, 1/beta) to create a 2 × 2 matrix
[,1] [,2]
[1,] 1 1
[2,] 1 1
rather than using c(beta[1]^2, 1/beta[2])) to create a vector c(1^2, 1/1) of length 2.
When you performed your matrix multiplication cbind(y, z) %*% ..., you consequently multiplied your 10 × 2 matrix cbind(y, z) by a 2 × 2 matrix, which yielded a 10 × 2 matrix as the output for your function u(). With a properly generated vector, however, the product would have been a 10 × 1 matrix.
Unsurprisingly, numDeriv::jacobian() gave a different result for a 10 × 2 matrix than for your intended 10 × 1 matrix.
I can give you a generalized function h(), which can be wrapped by u() to create the "pseudofunction" you describe here:
function(beta) c(y[1] * beta[1]^2 + z[1] * 1/beta[2], y[2] * beta[1]^2 + z[2] * 1/beta[2], : : : : y[i] * beta[1]^2 + z[i] * 1/beta[2])
For h(), we supply the arguments
beta: a numeric vector, of length n.funs: a list of n functions....: n numeric vectors each of length m, which will be consolidated as columns in a single m × n matrix A. Alternatively, the numeric matrix A itself.expand: a logical value indicating how funs should be applied to beta, to yield the structure by which the m × n matrix A will be multiplied:
TRUE: Apply to beta (as a whole) each of the n listed funs, and then consolidate each of the n results as a column of length n in an n × n matrix B.FALSE: Apply the the ith function in funs to the ith element in beta, and consolidate each of the n results as an element in a vector b of length n.and we receive either the m × n matrix AB (expand = TRUE) or the vector Ab of length m (expand = FALSE). Your purposes require the latter as input to pracma::jacobian().
Here is the definition of h()
h <- function(beta, funs, ..., expand = FALSE) {
# If there is only one function, encapsulate it in a list for mapply.
if(!is.list(funs)) {
funs <- list(funs)
}
# If expansion is desired, encapsulate beta in a list for mapply, to yield
# a set of vectors that can be consolidated as columns into a matrix.
# Otherwise, do neither, to yield a set of numbers consolidated as elements
# in a vector.
if(isTRUE(expand)) {
beta <- list(beta)
consolidate <- cbind
} else {
beta <- as.vector(beta)
consolidate <- base::c
}
return(
as.matrix(cbind(...) %*%
do.call(consolidate,
mapply(FUN = function(f, x) {
as.vector(sapply(X = x, FUN = f, simplify = TRUE))
},
funs, beta,
SIMPLIFY = FALSE)))
)
}
and here is the convenience function u() to wrap h() for your specific purposes:
y <- 1:10
z <- 10:19
u <- function(beta) {
h(beta = beta, funs = list(function(x){x^2}, function(x){1/x}), y, z, expand = FALSE)
}
You can now use
pracma::jacobian(u, c(1,1))
to obtain your desired output:
[,1] [,2]
[1,] 2 -10
[2,] 4 -11
[3,] 6 -12
[4,] 8 -13
[5,] 10 -14
[6,] 12 -15
[7,] 14 -16
[8,] 16 -17
[9,] 18 -18
[10,] 20 -19