问题描述
假设我想构造以下函数：
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
[,] -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 的雅可比矩阵？
jacobian 函数的结构如下
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 -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
解决方法
注意
你在一个特定的地方出错了：
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 矩阵。
不出所料，对于 10 × 2 矩阵，numDeriv::jacobian() 给出的结果与预期的 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个 nfunction个。 -
...：n 个长度为 m 的数值向量，它们将被合并为一个 m × n 矩阵 A。或者，数字matrixA 本身。 -
expand：一个逻辑值，指示如何将funs应用于beta，以产生 m × n 矩阵 A 将相乘：-
TRUE：应用于beta（作为一个整体）列出的每个nfuns，然后合并每个n 结果为 n × n 矩阵 B 中长度为 n 的列。 -
FALSE：将function中的 ithfuns应用到beta中的 ith 元素，并将每个 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))
},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))
获得你想要的输出：
[,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