Создание аргументов функции из именованного списка (с применением к 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)

Я не могу повторно использовать свою функцию логарифмического правдоподобия - мне придется переписать ее, чтобы иметь параметр бета3.Я хотел бы сделать что-то вроде:

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

, если , тогда я мог бы как-то настроить аргумент start в mle, чтобы учесть мою теперь функцию логарифмического правдоподобия для векторного ввода,или, исключая это, иметь функцию, которая создает функцию правдоподобия во время выполнения, скажем, путем создания именованного списка аргументов и последующего использования его для определения функции, например, что-то вроде этого:

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

Вы можете использовать функцию mle2 из пакета bbmle, которая позволяет передавать векторы в качестве параметров.Вот пример кода.

# 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

Answer 2

Возможно, будет проще использовать optim напрямую;это то, что mle использует в любом случае.

ll2 <- function(par, X, Y){
  beta <- matrix(c(par[-1]), ncol=1)
  -sum(log(dnorm(Y - X %*% beta, 0, par[1])))
}
getp <- function(X, sigma=1, beta=0.1) {
  p <- c(sigma, rep(beta, ncol(X)))
  names(p) <- c("sigma", paste("beta", 0:(ncol(X)-1), sep=""))
  p
}

set.seed(5)
n <- 100
df <- data.frame(x1 = runif(n), x2 = runif(n), y = runif(n))
Y <- df$y
X1 <- model.matrix(y ~ x1, data = df) 
X2 <- model.matrix(y ~ x1 + x2, data = df)
optim(getp(X1), ll2, X=X1, Y=Y)$par
optim(getp(X2), ll2, X=X2, Y=Y)$par

С выводом

> optim(getp(X1), ll2, X=X1, Y=Y)$par
      sigma       beta0       beta1 
 0.30506139  0.47607747 -0.04478441 
> optim(getp(X2), ll2, X=X2, Y=Y)$par
      sigma       beta0       beta1       beta2 
 0.30114079  0.39452726 -0.06418481  0.17950760 

Answer 3

Возможно, это не то, что вы ищете, но я бы сделал это следующим образом:

mle2(y ~ dnorm(mu, sigma),parameters=list(mu~x1 + x2), data = df,
    start = list(mu = 1,sigma = 1))

mle2(y ~ dnorm(mu,sigma), parameters = list(mu ~ x1), data = df,
    start = list(mu=1,sigma=1))

Возможно, вы сможете адаптировать эту формулировку для полинома, хотя dmultinom может не работать- вам может понадобиться написать Dmultinom(), который взял бы матрицу полиномиальных выборок и вернул (log) вероятность.

Answer 4

Код R, предоставленный Рамнатхом, также может быть применен к функции optim, поскольку также принимает векторы в качестве параметров.