r - 从命名列表创建函数参数(使用 stats4::mle 的应用程序)

Question

我应该首先说明我想要做什么:我想使用 mle 函数，而不必每次我想尝试不同的模型规范时都重新编写我的对数似然函数。因为 mle 需要一个指定的起始值列表，所以您显然不能只将对数似然函数编写为采用参数向量。一个简单的例子:

假设我想通过最大似然拟合线性回归模型，首先，我忽略了一个预测变量:

n <- 100
df <- data.frame(x1 = runif(n), x2 = runif(n), y = runif(n))
Y <- df$y
X <- model.matrix(lm(y ~ x1, data = df))

# define log-likelihood function 
ll <- function(beta0, beta1, sigma){
  beta = matrix(NA, nrow=2, ncol=1)
  beta[,1] = c(beta0, beta1)
  -sum(log(dnorm(Y - X %*% beta, 0, sigma)))
}
library(stats4)
mle(ll, start = list(beta0=.1, beta1=.2, sigma=1)

现在，如果我想适应不同的模型，请说:

m <- lm(y ~ x1 + x2, data = df)

我无法重复使用我的对数似然函数——我必须重新编写它才能具有 beta3 参数。我想做的是这样的:

ll.flex <- function(theta){
 # theta is a vector that I can use directly 
  ...
}

如果然后我可以以某种方式调整 mle 中的起始参数来解释我现在的向量输入对数似然函数，或者禁止这种情况，有一个在运行时构造对数似然函数的函数-time，例如通过构造命名参数列表，然后使用它来定义函数，例如，如下所示:

X <- model.matrix(lm(y ~ x1 + x2, data = df))
arguments <- rep(NA, dim(X)[2])
names(arguments) <- colnames(X)
ll.magic <- function(bring.this.to.life.as.function.arguments(arguments)){...} 

更新:

我最终编写了一个辅助函数，它可以将任意数量的命名参数 x1、x2、x3... 添加到传递的函数 f 中。

add.arguments <- function(f,n){
  # adds n arguments to a function f; returns that new function 
  t = paste("arg <- alist(",
  paste(sapply(1:n, function(i) paste("x",i, "=",sep="")), collapse=","),
  ")", sep="")
  formals(f) <- eval(parse(text=t))
  f
 }

它很丑陋，但它完成了工作，让我可以动态地重构我的对数似然函数。

Answer 1

您可以使用 bbmle 包中的 mle2 函数，它允许您将向量作为参数传递。这是一些示例代码。

# REDEFINE LOG LIKELIHOOD
ll2 = function(params){
  beta = matrix(NA, nrow = length(params) - 1, ncol = 1)
  beta[,1] = params[-length(params)] 
  sigma    = params[[length(params)]]
  minusll  = -sum(log(dnorm(Y - X %*% beta, 0, sigma)))
  return(minusll)
}

# REGRESS Y ON X1
X <- model.matrix(lm(y ~ x1, data = df))
mle2(ll2, start = c(beta0 = 0.1, beta1 = 0.2, sigma = 1), 
  vecpar = TRUE, parnames = c('beta0', 'beta1', 'sigma'))

# REGRESS Y ON X1 + X2

X <- model.matrix(lm(y ~ x1 + x2, data = df))
mle2(ll2, start = c(beta0 = 0.1, beta1 = 0.2, beta2 = 0.1, sigma = 1), 
      vecpar = TRUE, parnames = c('beta0', 'beta1', 'beta2', 'sigma'))

这给你

Call:
mle2(minuslogl = ll2, start = c(beta0 = 0.1, beta1 = 0.2, beta2 = 0.1, 
    sigma = 1), vecpar = TRUE, parnames = c("beta0", "beta1", 
    "beta2", "sigma"))

Coefficients:
     beta0      beta1      beta2      sigma 
 0.5526946 -0.2374106  0.1277266  0.2861055