Как передать Якобиан функции nlsLM() в R?

Содержание

Вопрос или проблема
Ответ или решение
Проблема и её решение
Корректное использование Якобиана
Объяснение кода

Вопрос или проблема

Я работаю с нелинейным моделированием, используя пакет minpack.lm, в частности функцию nlsLM(). Я хочу предоставить пользовательский Якобиан для улучшения производительности модели, но постоянно сталкиваюсь со следующей ошибкой при попытке его включить:

Ошибка в model.frame.default(formula = ~N + TT, jac = fcn.jac) :
недопустимый тип (closure) для переменной '(jac)'

Я использую nlsLM(), потому что мне нужно создать объект модели для дальнейшего использования при построении прогнозов. Ниже представлена немного измененная версия примера из документации.

Вопрос: Как правильно передать Якобиан в функцию nlsLM()? Что я делаю не так?

library(minpack.lm)

# Значения, по которым будет имитироваться данные
x <- seq(0,5, length = 100)

# Модель, основанная на списке параметров
getPred <- function(parS, xx) parS$a * exp(xx * parS$b) + parS$c

# Значения параметров, используемые для имитации данных
pp <- list(a = 9, b = -1, c = 6)

# Смоделированные данные с шумом
simDNoisy <- getPred(pp, x) + rnorm(length(x), sd = .1)

# Подгонка модели с использованием nlsLM
mod <- nlsLM(simDNoisy ~ a * exp(b * x) + c,
             start = c(a = 3, b = -0.001, c = 1),
             trace = TRUE)

# Построение графика данных
plot(x, simDNoisy, main = "данные")

# Построение графика подогнанных значений
lines(x, fitted(mod), col = 2, lwd = 2)

# Теперь давайте введем другой пример, используя пользовательский Якобиан

# Создание убывающей косинусоидальной функции с шумом
TT <- seq(0, 8, length = 501)
tau <- 2.2
N0 <- 1000
a <- 0.25
f0 <- 8
Ndet <- N0 * exp(-TT/tau) * (1 + a * cos(f0 * TT))
N <- Ndet + rnorm(length(Ndet), mean = Ndet, sd = .01 * max(Ndet))

# Функция для имитации данных
f <- function(TT, tau, N0, a, f0) {
  expr <- expression(N0 * exp(-TT/tau) * (1 + a * cos(f0 * TT)))
  eval(expr)
}

# Аналитический градиент (Якобиан)
j <- function(TT, tau, N0, a, f0) {
  expr <- expression(N0 * exp(-TT/tau) * (1 + a * cos(f0 * TT)))
  c(eval(D(expr, "tau")), eval(D(expr, "N0")),
    eval(D(expr, "a")), eval(D(expr, "f0")))
}

# Определение функции остатка
fcn <- function(p, TT, N, fcall, jcall) {
  N - do.call("fcall", c(list(TT = TT), as.list(p)))
}

# Определение аналитического Якобиана
fcn.jac <- function(p, TT, N, fcall, jcall) {
  -do.call("j", c(list(TT = TT), as.list(p)))
}

# Попытка подогнать модель с пользовательским Якобианом
mod <- nlsLM(N ~ N0 * exp(-TT/tau) * (1 + a * cos(f0 * TT)),
             start = c(tau = 2.2, N0 = 1500, a = 0.25, f0 = 10),
             trace = TRUE, jac = fcn.jac)

Когда бы я ни пытался передать Якобиан (fcn.jac) в nlsLM(), это приводит к ошибке. Буду признателен за любые советы по правильному способу передачи Якобиана в функцию nlsLM()!

Ответ или решение

Для корректной передачи Якобиана в функцию nlsLM() из пакета minpack.lm в R важно следовать определённым правилам. Давайте разберёмся с вашей ошибкой, а также посмотрим, как правильно прописать собственный Якобиан.

Проблема и её решение

Ошибка, с которой вы столкнулись:

Error in model.frame.default(formula = ~N + TT, jac = fcn.jac) :
invalid type (closure) for variable '(jac)'

указана на то, что аргумент jac не распознаётся как функция для передачи Якобиана.

Для передачи Якобиана в nlsLM() необходимо, чтобы функция-Якобиан принимала параметры в определённом формате. Она должна принимать в качестве аргументов вектор с параметрами (например, p), переменные, для которых вы хотите вычислить Якобиан, и дополнительные параметры, если это необходимо.

Корректное использование Якобиана

Вот пример исправленной реализации с вашим исходным кодом:

library(minpack.lm)

# Создаём данные
TT <- seq(0, 8, length = 501)
tau <- 2.2
N0 <- 1000
a <- 0.25
f0 <- 8
Ndet <- N0 * exp(-TT / tau) * (1 + a * cos(f0 * TT))
N <- Ndet + rnorm(length(Ndet), mean = 0, sd = 0.01 * max(Ndet))

# Функция для симуляции данных
f <- function(TT, tau, N0, a, f0) {
    N0 * exp(-TT / tau) * (1 + a * cos(f0 * TT))
}

# Функция для вычисления Якобиана
j <- function(TT, tau, N0, a, f0) {
    matrix(c(
      N0 * exp(-TT/tau) * (1 + a * cos(f0 * TT)) * TT / (tau^2),         # d/dtau
      exp(-TT/tau) * (1 + a * cos(f0 * TT)),                             # d/dN0
      N0 * exp(-TT/tau) * cos(f0 * TT) * TT,                             # d/da
      -N0 * exp(-TT/tau) * sin(f0 * TT) * TT                             # d/df0
    ), nrow = length(TT), ncol = 4, byrow = FALSE)
}

# Определяем функцию невязки
fcn <- function(p, TT, N) {
    N - f(TT, p[1], p[2], p[3], p[4])
}

# Определяем аналитику Якобиана
fcn.jac <- function(p, TT, N) {
    -j(TT, p[1], p[2], p[3], p[4])
}

# Попытка подгонки модели с кастомным Якобианом
mod <- nlsLM(N ~ f(TT, tau, N0, a, f0),
             start = c(tau = 2.2, N0 = 1500, a = 0.25, f0 = 10),
             trace = TRUE, jac = fcn.jac)

Объяснение кода

Функция f — определяет, как ваши параметры (tau, N0, a, f0) влияют на обрабатываемые данные.
Функция j — это якобиан, который вычисляет частные производные вашей функции по каждому из параметров.
Функция fcn — вычисляет невязку между наблюдаемыми и модельными значениями.
Функция fcn.jac — возвращает Якобиан, необходимые для оптимизации, как и было указано в nlsLM().

Таким образом, исправив реализацию вашего кода, вы сможете корректно передать Якобиан в функцию nlsLM(), что должно значительно улучшить производительность вашей модели.