Нелинейная оптимизация в R с Nloptr против Excel

Кажется, я не могу найти ответ для своей конкретной проблемы. Я пытаюсь минимизировать среднюю абсолютную ошибку между вектором действительного значения и линейной комбинацией моделей следующим образом:

library(nloptr)

df <- data.frame(
  real = c(24.2418, 21.7374, 7.203, 115.0233, 16.632, 5.4644, 27.8917, 0.5904, 0.633, 105.3656, 110.0564, 122.9399, 23.0418, 4.2186, 109.5453), 
  m_1 = c(17.7790979854662, 32.93639375, 6.94294375000066, 98.909065625, 80.1068562499999, 11.656556250002, 39.868921875, 0.859157480988586, 0.612625482376216, 112.580383416359, 151.896765176957, 155.81521460987, 7.3257, 4.1268, 41.5711703205879), 
  m_2 = c(25.5062900474607, 32.7709877317137, 6.94317649489279, 98.898593448344, 80.1114267521597, 11.6563450007055, 39.8691722409057, 0.907260912997819, 0.602795676399197, 114.183526809465, 139.51052151724, 149.993624420536, 6.85002142907728, 4.66668862149305, 70.7527906311631), 
  m_3 = c(27.1495912199999, 40.2339353399999, 7.10339542, 87.1444967133334, 58.4563384933332, 11.1378198366666, 37.6030141333333, 0.852288459999999, 0.681724560000002, 100.101136186666, 118.536525109999, 136.923264319999, 5.64763034333333, 4.8659515, 70.12675835), 
  m_4 = c(25.511590625, 32.9363937499999, 7.00050769504031, 98.3660974929738, 80.10685625, 11.65655625, 39.868921875, 0.665580984791989, 0.498756215272925, 85.791042265746, 135.619508469251, 140.946144588455, 5.05824305930683, 3.25333636845094, 22.0908752674237), 
  m_5 = c(25.6118152340655, 34.5646719234769, 6.82416840383483, 91.5582383465651, 84.4603847826215, 11.3405620277701, 40.7906062254877, 0.908706704665592, 0.602817399156822, 114.326905157898, 139.595783699511, 150.046375909198, 6.8511793011574, 4.6622942290559, 56.2753212961812), 
  m_6 = c(21.9868574376585, 44.3147731773466, 6.38315486686481, 100.303757097094, 9.13921010739697, 7.83817900918309, 31.5458855316741, 1.09960505333834, 0.817751834425678, 101.110814846224, 145.55847538105, 142.82362305075, 7.61732986965459, 4.6774198307473, 67.5821464371521)
)

best_dist <- function(x) {
  output <- df$m_1 * x[1] + df$m_2 * x[2] + df$m_3 * x[3] + 
    df$m_4 * x[4] + df$m_5 * x[5] + df$m_6 * x[6]
  mean(abs(output - df$real))
}

restriction <- function(x) sum(x) - 1

nloptr(
  x0 = rep(1 / 6, 6), 
  eval_f = best_dist, 
  lb = rep(0, 6), 
  ub = rep(1, 6), 
  eval_g_eq = restriction, 
  opts = list(algorithm = "NLOPT_GN_ISRES", xtol_rel = 1e-16, maxeval = 1e4)
)

Как вы могли прочитать, я использую nloptr пакет. Приведенный выше код дает неоптимальный результат 14,85 для целевой функции, и все параметры являются начальными параметрами. Вы можете изменить начальные параметры на какой-то другой вектор и все равно не получите оптимального решения.

Однако с помощью решателя Excel можно легко получить результат 10,77 для целевой функции и (0, 0, .15, 0, 0, .85) для параметров.

Я попытался использовать алгоритм с градиентом, но я не могу понять синтаксис правильно. Вот моя другая попытка.

gradient <- function(x) {
  output <- df$m_1 * x[1] + df$m_2 * x[2] + df$m_3 * x[3] + 
    df$m_4 * x[4] + df$m_5 * x[5] + df$m_6 * x[6]
  err <- output - df$real

  c(
    - sum(sign(err) * df$m_1), 
    - sum(sign(err) * df$m_2), 
    - sum(sign(err) * df$m_3), 
    - sum(sign(err) * df$m_4), 
    - sum(sign(err) * df$m_5), 
    - sum(sign(err) * df$m_6)
  )
}

nloptr(
  x0 = runif(6), 
  eval_f = best_dist, 
  eval_grad_f = gradient, 
  lb = rep(0, 6), 
  ub = rep(1, 6), 
  eval_g_eq = restriction, 
  opts = list(algorithm = "NLOPT_GN_ISRES", xtol_rel = 1e-16, maxeval = 1e4)
)
Источник

1 ответ

Это похоже на регрессию LAD с боковым ограничением. LP (линейное программирование) формулировки для регрессии LAD показаны здесь. Используя ваше наименование данных (real[i], m[i,j]), мы имеем:

min sum(i, | real[i] - sum(j,m[i,j]*x[j]) |  ) 
subject to
   sum(j, x[j]) = 1
   0 <= x[j] <= 1

Это может быть линеаризовано (с использованием разделения переменных) следующим образом:

min sum(i, r1[i]+r2[i])
subject to
   r1[i]-r2[i] = real[i] - sum(j,m[i,j]*x[j])
   sum(j, x[j]) = 1
bounds:
   r1[i], r2[i] >= 0       (positive variables)
   0 <= x[j] <= 1

Это чистый LP и может быть решен с помощью любого решателя LP.

Источник
Другие вопросы по тегам