Оптимизация с несколькими параметрами в R

ProgramBox

Оптимизация с несколькими параметрами в R

Post author:admin
Запись опубликована:3 апреля, 2023
Post category:Вопросы по программированию

Вопрос:

У меня есть несколько точек данных для каждого месяца:

 datapoints = c(0.166247133507398, 0.425481263534677, 0.411408800035613, 0.435293979295827, 
0.423245332292126, 0.346365073582275, 0.425199833179369, 0.432137205696319, 
0.43081430546139, 0.432920308468688, 0.453136718369095, 0.455262858225237, 
0.455296396323467, 0.463279988479843, 0.476385172463429, 0.249218233044765, 
0.471649455997085, 0.478360455837937, 0.476885553827009, 0.478454365885292, 
0.469363359162398, 0.482051012672114, 0.472681867759822, 0.473375139224335, 
0.466772123979772, 0.470227689772718, 0.31968277436218, 0.469224848893103, 
0.462242410659246, 0.460592456233742, 0.47620986576363, 0.479597725173361, 
0.463367143161687, 0.445818446073818, 0.448169906772321, 0.469654391229704, 
0.448246061325598, 0.371469823189851, 0.444850426677065, 0.458648112858549, 
0.447860339537274, 0.459419695588122, 0.460269819008956, 0.47749167939716, 
0.473072390305622, 0.476016821990266, 0.4568557470378, 0.449672268375193, 
0.380564883000597, 0.444441519806569, 0.460583921380099, 0.461494033891503, 
0.469138527283594, 0.45150277521197, 0.46718546076104, 0.477345167601084, 
0.472688714814486, 0.487739306383874, 0.480842942430047, 0.391863568128969, 
0.480438541452507, 0.489438220266693, 0.476177921939346, 0.473645046232311, 
0.465869413362134, 0.455922096430803, 0.461787856186048, 0.467509464988477, 
0.484448645266681, 0.471468467400354, 0.395605537227465, 0.479981910572418, 
0.48446619841535, 0.435164680279421, 0.450079104167341, 0.44282935863006, 
0.459921867384553, 0.432157339753678, 0.457218806281876, 0.50047675208151, 
0.436606373021543, 0.40761740347071, 0.415360136419418, 0.427831358793833, 
0.432500363208447, 0.462355479936914, 0.427092963784101, 0.456745139905428, 
0.457223277757524, 0.459387550114517, 0.490044004170164, 0.436484845722895, 
0.420663824652525, 0.435432613404983, 0.432974547875276, 0.457827421432496, 
0.488379378067953, 0.482342084402802, 0.475817074700216, 0.432648247480694, 
0.396584389592281, 0.449595227767319, 0.457362567053931, 0)

 months= c(-1, -2, -3, -4, -5, -6, -7, -8, -9, -10, -11, -12, -13, -14, 
-15, -16, -17, -18, -19, -20, -21, -22, -23, -24, -25, -26, -27, 
-28, -29, -30, -31, -32, -33, -34, -35, -36, -37, -38, -39, -40, 
-41, -42, -43, -44, -45, -46, -47, -48, -49, -50, -51, -52, -53, 
-54, -55, -56, -57, -58, -59, -60, -61, -62, -63, -64, -65, -66, 
-67, -68, -69, -70, -71, -72, -73, -74, -75, -76, -77, -78, -79, 
-80, -81, -82, -83, -84, -85, -86, -87, -88, -89, -90, -91, -92, 
-93, -94, -95, -96, -97, -98, -99, -100, -101, -102, -103)

Я хотел бы применить следующую функцию к точкам данных для оценки параметров a, b и c, которая стремится минимизировать сумму квадратов различий между исходными точками данных (кривой) и подобранной.

 Function -> a*exp((b c*(-numberofmonths ))/(-numberofmonths))

Я попытался использовать optimix в R следующим образом:

 library(optimx)

#Create function 
fittedCCF.f<- function(a,b,c){
  optimum<-a*exp((b c*(-months))/(-months))
  return(optimum)
}

#Run optimisation on function to determine a,b and c parameters
result=optimx(datapoints,fittedCCF.f)

Следующая ошибка показывает:

 Error in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower,  : 
  Cannot evaluate function at initial parameters

1. Можете ли вы опубликовать образцы данных? Пожалуйста, отредактируйте вопрос с выводом dput(datapoints) .

2. Правильно ли аналитическое выражение функции? Правильны ли скобки ? Обратите внимание, что в функции отсутствует независимая переменная, у нее есть только параметры. И, пожалуйста, включите вызов library для загрузки пакета, который optimx можно найти, поскольку это не базовая функция R.

3. Скобки в функции указаны правильно. Независимой переменной является количество месяцев. Выполняется вызов library(optimx)

Ответ №1:

Хотя уже есть принятый ответ, вот еще один с optimx . Может быть полезно иметь более общее решение, с несколькими методами, используемыми для нахождения минимума функции, и с двумя разными начальными точками, одна из которых является вектором, близким к нулю ( .Machine$double.eps^0.5 ), а другая — вектором 1 ‘s.

Установленная функция та же, переписанная, возможно, для большей читаемости. Вспомогательная функция convergence возвращает TRUE значение, convcode равное нулю, и если конечный градиент близок к нулю ( kkt1 равно TRUE ).

 library(optimx)

MSE.f <- function(x, target, months){
  a <- x[1]
  b <- x[2]
  c <- x[3]
  y <- a*exp(c - b/months)
  sum((target - y)^2)/length(target)
}

convergence <- function(ans){
  ans[['convcode']] == 0 amp; ans[['kkt1']] amp; !is.na(ans[['kkt1']])
}

Это похоже на принятый ответ. В n приведенном ниже коде нет необходимости.

 length(datapoints)
# 104
length(months)
# 103
datapoints <- datapoints[seq_along(months)]

Теперь код оптимизации.

 p0 <- rep(.Machine$double.eps^0.5, 3)
o0 <- optimx(par = p0, MSE.f, 
             control = list(all.methods = TRUE),
             target = datapoints, months = months)

p1 <- c(1, 1, 1)
o1 <- optimx(par = p1, MSE.f, 
             control = list(all.methods = TRUE),
             target = datapoints, months = months)

i0 <- convergence(o0)
i1 <- convergence(o1)

o0[i0, ]
#                p1         p2           p3       value fevals gevals niter convcode kkt1  kkt2 xtime
#BFGS     0.4151490 -0.6022883  0.103498081 0.001354874     32     23    NA        0 TRUE FALSE 0.016
#L-BFGS-B 0.4512921 -0.6021743  0.020016393 0.001354874     13     13    NA        0 TRUE FALSE 0.003
#nlm      0.4585223 -0.6021536  0.004120473 0.001354874     NA     NA     7        0 TRUE FALSE 0.001
#nlminb   0.4572070 -0.6021818  0.006995190 0.001354874     14     40    10        0 TRUE FALSE 0.002
#spg      0.4704948 -0.6005138 -0.021720301 0.001354877     51     NA    43        0 TRUE FALSE 0.627
#ucminf   0.4572098 -0.6021802  0.006988472 0.001354874     10     10    NA        0 TRUE FALSE 0.001
#Rcgmin   0.4451012 -0.6021822  0.033829649 0.001354874     53     22    NA        0 TRUE FALSE 0.004
#Rvmmin   0.4129271 -0.6021820  0.108860350 0.001354874     31     20    NA        0 TRUE FALSE 0.007
#newuoa   0.4512979 -0.6021820  0.020003758 0.001354874    260     NA    NA        0 TRUE FALSE 0.005
#bobyqa   0.4812281 -0.6021820 -0.044210007 0.001354874    267     NA    NA        0 TRUE FALSE 0.006

o1[i1, ]
#                  p1         p2           p3       value fevals gevals niter convcode kkt1  kkt2 xtime
#BFGS     -13.5823818 -0.7876065 -13.58242467 0.202300550      4      2    NA        0 TRUE FALSE 0.000
#L-BFGS-B   0.2858336 -0.6021726   0.47672093 0.001354874     17     17    NA        0 TRUE FALSE 0.003
#nlm      -13.1595761 -1.0101522 -19.32529032 0.202285737     NA     NA    10        0 TRUE FALSE 0.001
#nlminb     0.2933911 -0.6021820   0.45062493 0.001354874     17     61    15        0 TRUE FALSE 0.002
#spg        0.4505517 -0.6006739   0.02159383 0.001354877     69     NA    59        0 TRUE FALSE 0.657
#ucminf     0.2933888 -0.6021851   0.45063208 0.001354874     15     15    NA        0 TRUE FALSE 0.001
#Rcgmin   -16.4515799 -1.1736319 -17.27751080 0.202286132    225    252    NA        0 TRUE FALSE 0.029
#newuoa     0.1723519 -0.6021814   0.98259278 0.001354874    467     NA    NA        0 TRUE FALSE 0.010
#bobyqa     0.5608663 -0.6021823  -0.19735125 0.001354874    540     NA    NA        0 TRUE FALSE 0.011

Похоже, что у начального решения, близкого к нулю, меньше проблем с конвергенцией, и решения, найденные несколькими методами, очень похожи. Начальная точка c(1, 1, 1) , напротив, показывает больший диапазон конечных значений. Значения функции также менее стабильны.

Ответ №2:

Прежде всего, ваши datapoints и months векторы должны быть одинаковой длины.

 length(datapoints)
# 104
length(months)
# 103

Итак, чтобы сделать вычисления возможными, я отбросил последнее значение datapoints .

 n <- length(datapoints)
datapoints <- datapoints[-n]

Тогда я не использовал optimx библиотеку, вот пример с базовой функцией R optim .

 opt_result <- optim(
    par = c(1, 1, 1)
    , fn = function(p) {
        a <- p[1]; b <- p[2]; c <- p[3]
        sq_err <- sum((datapoints - a*exp((b c*(-months))/(-months)))^2)/(n-1)
    }
    , method = "L-BFGS-B"
)

opt_result
# $par
# [1]  0.2858336 -0.6021726  0.4767209
# 
# $value
# [1] 0.001354874
# 
# $counts
# function gradient 
# 17       17 
# 
# $convergence
# [1] 0
# 
# $message
# [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH"

1. Большое вам спасибо! Кажется, что теперь это ближе к результату, однако тогда method = «L-BFGS-B», похоже, работает некорректно. Появляется следующая ошибка: Error in optim(par = c(1, 1, 1), fn = function(p) { : L-BFGS-B needs finite values of 'fn'

2. вероятно, это из-за большого значения внутри exp — если оно больше, log(.Machine$double.xmax) = 709.7827 значение exp станет Inf . Вы можете попробовать переписать свою формулу как log(datapoints) ~ log(a) log((b c*(-months))/(-months))

3. Было бы более уместно использовать вместо этого функцию nls.lm? rdocumentation.org/packages/minpack.lm/versions/1.2-1/topics /…

4. Согласно моему опыту работы с nls.lm , иногда это приводило к ошибкам, которые не могут сходиться, в других случаях работало довольно медленнее, чем optim с 'L-BFGS-B' методом, поэтому я не использовал его в производстве. Это было несколько лет назад, может быть, сегодня эта функция работает лучше 🙂