在函数中循环

问题描述

Jacobian 矩阵的封装

install.packages('pracma')
library('pracma')

假设我有以下一组非线性函数

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

现在假设 i 变得非常大，写出 i 的所有函数非常耗时。

我怎样才能写出这样的函数？ 我尝试了以下方法：

f <- function(beta) cbind(y,z) %*% (c(beta^2,(beta^-1)),ncol = 1,nrow = 2)

当我应用这个函数来构造雅可比矩阵时，它不是我想要的格式

y = 1:10
z = 10:19
f <- function(beta) cbind(y,z) %*% matrix(cbind(beta^2,1/beta),nrow = 2)
jacobian(f,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

再试一次：

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

给出以下结果：

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

第一列的前 10 个观察值和第二列的 11:20 观察值是我需要的一次。

有谁知道如何获得这两个向量？我应该在函数内部创建一个循环吗？

解决方法

您可以按如下方式在 R 中使用 vectorisation：

y <- 1:10
z <- 10:19
beta <- c(1,1)

f <- function(y,z,beta = c(1,1)) {
  y * beta[[1]] ** 2 + z * 1 / beta[[2]]
}

f(y,beta)

#> f(y,beta)
# [1] 11 13 15 17 19 21 23 25 27 29

编辑：

感谢 @Hans W. 的评论和对我的回答的改进。这是编辑：

# define y & z in global working environment:
y <- 1:10
z <- 10:19

# The function:
f <- function(beta) {
  y * beta[[1]] ** 2 + z * 1 / beta[[2]]
}

# Find Jacobian:
pracma::jacobian(f,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