时间戳

首页 > 解决方案 > R中n个变量的m个函数

r - R中n个变量的m个函数

问题描述

假设我要构造以下函数：

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

标签： rfunctionmatrixnonlinear-functions

解决方案

笔记

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

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]))用来创建c(1^2, 1/1)长度为 2 的向量。

当您执行矩阵乘法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: 一个listn秒 function。
  • ...: n 个长度为m的数字向量，它们将合并为单个m × n矩阵A中的列。或者，数字matrix A本身。
  • expand: 一个逻辑值，指示funs应如何应用于beta，以产生m × n矩阵A将乘以的结构：
    • TRUE：应用于beta（作为一个整体）列出的n 个funs中的每一个，然后将n 个结果中的每一个合并为n × n矩阵B中长度为n的列。
    • FALSE：将第i个functioninfuns应用于 中的第i个元素beta，并将n 个结果中的每一个合并为长度为n的向量b中的一个元素。

我们收到m × n矩阵AB ( ) 或长度为m ( )expand = TRUE的向量Ab。您的目的需要后者作为.expand = FALSEpracma::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

推荐阅读