使用 apply 系列函数编写对数似然。业绩损失？

Question

我正在使用 R 中的 apply() 系列函数，并尝试使用 apply() 编写一个对数似然函数。

这是假设高斯干扰的线性回归模型的对数似然：

# Likelihood function for the standard linear regression model

logL <- function(theta, data){
    # Return minus the log likelihood function for the standard linear regression model
    # y: endogenous variable
    # x: matrix of regressors
    y <- data[, 1]
    x <- data[, -1]
    N <- nrow(data)
    # This is the contribution to the log-likelihood of individual i. Initialized at 0.
    contrib <- 0
    beta <- head(theta, -1) # Every element but the last one
    sigma <- tail(theta, 1) # Only the last element
    for (i in 1:N){
        contrib <- contrib + (y[i] - beta%*%x[i,])**2
    }
    sigma <- abs(sigma)
    L <- -(1/(2*sigma^2)*contrib) - 1/2 * N * log(2*pi) - N * log(sigma)
    return(-L)
}

下面我们模拟一些数据，最小化负对数似然（相当于最大化对数似然）。

# Simulate some data
N <- 1000
x <- cbind(1, rnorm(N,0,sd=1), rnorm(N, 0, sd=2)) 

true_theta <- c(2, 3, 2, 4)

y <- true_theta[1:3]%*%t(x) + rnorm(N, mean = 0, sd = true_theta[4])

my_data <- cbind(t(y),x)

optim(c(1,1,1, 1), fn = logL, data = my_data, 
method = "L-BFGS-B",upper = c(Inf, Inf, Inf), lower=c(-Inf, -Inf, 0.01))

到目前为止一切顺利，我们得到的结果与用于模拟数据的结果相同。通过使用 rbenchmark 包，我得到优化步骤的 10 次复制在我的计算机上大约需要 4 秒。

benchmark(optim(c(1,1,1, 1), fn = logL, data = my_data,
method = "L-BFGS-B",upper = c(Inf, Inf, Inf), lower=c(-Inf, -Inf, 0.01)),
replications=10)

现在我尝试用 apply 函数替换 for 循环。为此，我将 contrib 定义为一个函数：

contrib <- function(beta, one_obs){
    y <- one_obs[1]
    x <- one_obs[-1]
    return((y - beta%*%x)**2)
}

还有新的对数似然函数：

logL2 <- function(theta, data){
# Return minus the log likelihood function for the standard linear   regression model
# y: endogenous variable
# x: matrix of regressors
N <- nrow(data)
beta <- head(theta, -1) # Every element but the last one
sigma <- tail(theta, 1) # Only the last element
sigma <- abs(sigma)
L <- -(1/(2*sigma^2)*sum(apply(data, FUN=contrib, beta = beta, 1))) 
- 1/2 * N * log(2*pi) - N * log(sigma)
return(-L)
}

这几乎是原来的两倍。现在，我可能误解了 apply 系列函数的作用，因为它们应该用于代码清晰而不是性能。但是，它们不应该比 for 循环慢，对吧？那么我的代码发生了什么？是否正在进行某种类型转换？我检查并 logL 返回一个矩阵， logL2 返回一个数字。我尝试使用 vapply() 因为它允许指定返回对象的类型，但 vapply() 似乎通过将每一列堆叠在一起将我的数据矩阵转换为向量。这会导致 contrib 函数不再工作：

logL2 <- function(theta, data){
# Return minus the log likelihood function for the standard linear regression model
# y: endogenous variable
# x: matrix of regressors
N <- nrow(data)
beta <- head(theta, -1) # Every element but the last one
sigma <- tail(theta, 1) # Only the last element
sigma <- abs(sigma)
L <- -(1/(2*sigma^2)*sum(vapply(data, FUN=contrib, beta = beta, FUN.VALUE = matrix(1)))) - 1/2 * N * log(2*pi) - N * log(sigma)
return(-L)
}

这就是我得到的结果：

class(logL2(theta = c(1,2,2,2), my_data))

Error in beta %*% x : non-conformable arguments 

那么我如何使用 apply 系列函数来使我的代码更具可读性，并且至少与使用 for 循环一样快？

Answer 1

您可以通过考虑 for 循环中涉及的数学来简化代码。

你的 for 循环是

contrib <- contrib + (y[i] - beta%*%x[i,])**2

现在这与计算所有(y[i] - beta %*% x[i, ])^2 并将它们相加相同。考虑beta %*% x[i, ]，您正在对 1x3 矩阵 (beta) 与 3x1 (x[i, ]) 进行矩阵乘法，得到 1x1 结果。所以你正在做的是将beta 矩阵乘以x 的每一行。 但是，使用矩阵乘法，您无论如何都可以同时进行所有操作，并得到一个 Nx1 矩阵！

即beta (1x3) %*% x (3xN) 会给你一个 1xN 的矩阵，然后从 y 中减去这个矩阵，y 也是一个长度为 N 的向量，独立地对每个差值求平方并将它们相加。这相当于你的 for 循环。

唯一的问题是你的 x 是 Nx3 而不是 3xN，所以我们先 t() 它：

contrib <- sum((y - beta %*% t(x))^2)

这完全消除了你的 for 循环。

logL2 <- function(theta, data){
    y <- data[, 1]
    x <- data[, -1]
    N <- nrow(data)
    beta <- head(theta, -1) # Every element but the last one
    sigma <- tail(theta, 1) # Only the last element
    contrib <- sum((y - beta %*% t(x))^2)
    sigma <- abs(sigma)
    L <- -(1/(2*sigma^2)*contrib) - 1/2 * N * log(2*pi) - N * log(sigma)
    return(-L)
}

library(rbenchmark)
benchmark(
    orig={orig.answer <- optim(c(1,1,1, 1), fn = logL, data = my_data,
method = "L-BFGS-B",upper = c(Inf, Inf, Inf), lower=c(-Inf, -Inf, 0.01))},
    new={new.answer <- optim(c(1,1,1, 1), fn = logL2, data = my_data,
method = "L-BFGS-B",upper = c(Inf, Inf, Inf), lower=c(-Inf, -Inf, 0.01))},
replications=10
)

产生

  test replications elapsed relative user.self sys.self user.child sys.child
2  new           10   0.306     1.00     0.332    0.048          0         0
1 orig           10   4.584    14.98     4.588    0.000          0         0

让我们检查一下我们没有犯错

all.equal(orig.answer, new.answer)
# [1] TRUE

作为一个风格点，为什么不将y 作为logL2 的第三个参数（而不是在开始时将cbind 设置为data，然后必须选择适当的行/列所有时间）？这使您不必一直执行y <- data[, 1] 和x <- data[, -1]。 IE。执行logL <- function (theta, x, y) { ... } 之类的操作，然后在您的optim() 调用中，您可以提供x 和y 参数而不是my_data。您甚至可以通过在开始时执行 t(x) 来获得进一步的改进（例如，在您调用 optim 时）因此不必每次调用 logL2 时都执行此操作？

logL3 <- function(theta, x, y){
  N <- length(y)
  beta <- head(theta, -1) # Every element but the last one
  sigma <- tail(theta, 1) # Only the last element
  contrib <- sum((y - beta %*% x)^2)
  sigma <- abs(sigma)
  L <- -(1/(2*sigma^2)*contrib) - 1/2 * N * log(2*pi) - N * log(sigma)
  return(-L)
}

benchmark(
  new=optim(c(1,1,1, 1), fn = logL2, data = my_data,
            method = "L-BFGS-B",upper = c(Inf, Inf, Inf), lower=c(-Inf, -Inf, 0.01)),
  new.new=optim(c(1,1,1, 1), fn = logL3, x=t(x), y=y,
            method = "L-BFGS-B",upper = c(Inf, Inf, Inf), lower=c(-Inf, -Inf, 0.01)),
  replications=100
)
     test replications elapsed relative user.self sys.self user.child sys.child
1     new          100   3.149    2.006     3.317    0.700          0         0
2 new.new          100   1.570    1.000     1.488    0.344          0         0

它的速度大约是原来的两倍。一般来说，如果你可以做一次而不是每次调用logL2（例如t(x)、data[, 1] 等），它会为你节省一些时间。

---

然而，关于您的原始问题（特别是与 *apply 函数有关：

• vapply 将 list 作为输入，而您的data 是一个矩阵，因此contrib 一次对data 的一个元素进行操作。 IE。 contrib 将 x 视为单个数字。因此不符合矩阵，因为您的矩阵乘法是将beta（1x3）与x（1x1）相乘，并且要使矩阵乘法起作用，您需要beta的列数等于@的行数987654366@。要使用vapply，您需要类似

  vapply(1:nrow(data), function(i) contrib(beta, data[i, ]), FUN.VALUE=1)
  

• （！我没有通过基准测试或其他任何方法测试这些语句。这正是我在经验中发现的）：在所有*apply 函数中，我发现apply() 很慢（通常比for 慢-环形）。它对于代码的简洁性很方便（“为每一行执行此操作”，或“为每一列执行此操作”类型的任务：而不是很多 data[i, ] 它只是 apply(.., MARGIN=1)），但是如果您需要速度，请执行循环或使用其他表亲之一，例如 vapply、lapply 或 sapply。
• vapply、lapply 很快。 sapply 也是如此，但通常前两者之一更快（sapply 更易于使用，因为 vapply 的 FUN.VALUE 位正在为您制定。或者如果您知道 FUN.VALUE不会总是一样的，它等同于lapply，所以你不妨使用它。由于sapply 为你完成了所有这些工作，它可以更容易使用，但速度会稍慢）。
• 最快的是如果你可以使用一些数学来避免循环！例如如果你可以像我在这里所做的那样用矩阵乘法来改写你的循环。虽然这只适用于极少数情况。