R中n个变量的m个函数

Question

假设我要构造如下函数：

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])

假设我有以下数据。

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.

分别为 f 和 g 产生所需的输出：

     [,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

现在我可以合并这两个矩阵来获得 f 和 g 的雅可比。但是，我只想要一个函数来获得所需的输出。

我尝试了以下方法，但这并没有产生我想要的结果：

u <- function(beta) (cbind(y, z) %*% cbind(beta^2,1/beta))
jacobian(u, c(1,1))

给出不正确的输出：

     [,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

有谁知道如何组合函数 f 和 g 以获得 10 x 2 雅可比矩阵？

雅可比函数的结构如下

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.

我想要得到的输出是

     [,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

Answer 1

注意

你在一个特定的地方出错了：

u <- function(beta) (cbind(y, z) %*% cbind(beta^2,1/beta))
#                                    ^^^^^^^^^^^^^^^^^^^^
#                                            HERE

您使用cbind(beta^2, 1/beta) 创建了一个 2 × 2 矩阵

     [,1] [,2]
[1,]    1    1
[2,]    1    1

而不是使用c(beta[1]^2, 1/beta[2])) 创建一个长度为2 的向量c(1^2, 1/1)。

当您执行矩阵乘法 cbind(y, z) %*% ... 时，您将 10 × 2 矩阵 cbind(y, z) 乘以一个 2 × 2 矩阵，得到一个 >10 × 2 矩阵作为函数u() 的输出。然而，使用正确生成的向量，乘积将是一个 10 × 1 矩阵。

不出所料，numDeriv::jacobian() 为 10 × 2 矩阵提供的结果与您预期的 10 × 1 矩阵不同。

广义解

我可以给你一个泛化函数h()，它可以被u() 包装，以创建你在这里描述的“伪函数”：

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])

对于h()，我们提供参数

• beta: 一个数字向量，长度为 n。
• funs：list of n functions.
• ...: n 个数字向量，每个长度为 m，将合并为单个 m × n 中的列em> 矩阵A。或者，数字 matrix A 本身。
• expand：一个逻辑值，指示如何将funs 应用于beta，以产生 m × n 矩阵 A 将相乘：
  • TRUE：应用于beta（作为一个整体）列出的funs中的每一个n，然后将每一个n结果合并为一列长度n 在 n × n 矩阵 B 中。
  • FALSE：将funs中的第i个function应用到beta中的第i个元素，并合并每个n 结果是长度为 n 的向量 b 中的一个元素。

我们收到 m × n 矩阵 AB (expand = TRUE) 或向量 Ab长度 m (expand = FALSE)。您的目的需要后者作为pracma::jacobian() 的输入。

这里是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)))
  )
}

这里是方便函数 u() 为您的特定目的包装 h()：

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)
}

您现在可以使用

pracma::jacobian(u, c(1,1))

获得你想要的输出：

      [,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