R中n个变量的m个函数
m functions of n variables in R
假设我要构造以下函数：
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
备注
你错在一处：
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: n functions.list....：n 个长度为 m 的数值向量，它们将合并为单个 [= 中的列79=]m×n矩阵A。或者，数字 matrix A 本身。expand：一个逻辑值，指示如何将 funs 应用于 beta，以产生 m × n 矩阵 A 将相乘：
TRUE：应用于beta（作为一个整体）列出的n中的每一个funs，然后合并每一个n 结果作为长度为 n 的列在 n × n矩阵B.FALSE：将funs中的第i第function应用到i beta 中的第 th 个元素，并将每个 n 结果合并为长度为 的向量 b 中的一个元素n.
并且我们收到 m × n 矩阵 AB (expand = TRUE) 或长度为 m (expand = FALSE) 的向量 Ab。您的目的需要后者作为 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
假设我要构造以下函数：
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
备注
你错在一处：
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: nfunctions.list....：n 个长度为 m 的数值向量，它们将合并为单个 [= 中的列79=]m×n矩阵A。或者，数字matrixA 本身。expand：一个逻辑值，指示如何将funs应用于beta，以产生 m × n 矩阵 A 将相乘：TRUE：应用于beta（作为一个整体）列出的n中的每一个funs，然后合并每一个n 结果作为长度为 n 的列在 n × n矩阵B.FALSE：将funs中的第i第function应用到ibeta中的第 th 个元素，并将每个 n 结果合并为长度为 的向量 b 中的一个元素n.
并且我们收到 m × n 矩阵 AB (expand = TRUE) 或长度为 m (expand = FALSE) 的向量 Ab。您的目的需要后者作为 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