J'ai un certain nombre de points de données pour chaque mois:

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)

Je voudrais ajuster la fonction suivante sur les points de données pour estimer les paramètres a, b et c qui cherchent à minimiser la somme des différences carrées entre les points de données d'origine (courbe) et celui ajusté.

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

J'ai essayé d'utiliser optimix dans R comme suit:

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)

L'erreur suivante s'affiche:

Error in optimx.check(par, optcfg$ufn, optcfg$ugr, optcfg$uhess, lower,  : 
  Cannot evaluate function at initial parameters
r
1
Monica Odysseos CY 2 sept. 2020 à 14:46

2 réponses

Meilleure réponse

Bien qu'il y ait déjà une réponse acceptée, en voici une autre avec optimx. Il peut être utile d'avoir une solution plus générale, avec plusieurs méthodes utilisées pour trouver le minimum de la fonction et avec deux points de départ différents, l'un un vecteur proche de zéro (.Machine$double.eps^0.5) et l'autre un vecteur de {{X2} }.

La fonction ajustée est la même, réécrite pour être peut-être plus lisible. Une fonction auxiliaire convergence renvoie TRUE si le convcode est nul et si le dégradé final est proche de zéro (kkt1 est 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 & ans[['kkt1']] & !is.na(ans[['kkt1']])
}

C'est comme dans la réponse acceptée. Il n'y a pas besoin de n dans le code qui suit.

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

Maintenant, le code d'optimisation.

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

Une solution initiale proche de zéro semble avoir moins de problèmes de convergence et les solutions trouvées par les différentes méthodes sont très similaires. Le point initial c(1, 1, 1), au contraire, montre une plage plus large dans les valeurs finales. Les valeurs de la fonction sont également moins stables.

2
Rui Barradas 2 sept. 2020 à 13:56

Tout d'abord, vos vecteurs datapoints et months doivent être de la même longueur.

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

Donc, pour rendre le calcul possible, j'ai écarté la dernière valeur de datapoints.

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

Ensuite, je n'ai pas utilisé la bibliothèque optimx, voici un exemple avec la fonction de base 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"
3
inscaven 2 sept. 2020 à 12:27