Bonjour à tous,
Débutant en programmation sous R je me retrouve bloqué pour l'optimisation de paramètres avec optim() mais aussi avec l'utilisation de ode() du package DeSolve.
Explication:
Je cherche à optimiser une courbe de croissance linéaire mais avec 2 phases (donc 2 taux de croissance : a et b) et une valeur de transition Mc
j'ai donc besoin d'optimiser 3 paramètres : Mc (valeur transitoire), a (taux de croissance avant Mc) et b (taux de croissance après Mc), à l'aide d'un jeu de données. Voici mon code :
#library deSolve
library(deSolve)
#==============================================================
# valeurs initiales des paramètres Mc, a et b dans le modéle
#==============================================================
val.ini=c(1, 0.01,0.01)
#=========================
données utilisables pour essai
#=========================
yobs= c(2.42,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,3.39,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,4.28,
NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,5.19,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,6.96,
NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,8.48,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA)
#je dispose de données que tous les 14 jours, mais je fais tourner le model avec un pas de temps de 1 jour, d'où les NA
#==============================================================
# fonction à optimiser
#==============================================================
croissance=function(t,y,parms){
with(as.list(c(y,parms)),{
if(y<Mc)
{dy=a}
else {dy=b}
list(c(dy))
}) }
# fonction d'optimisation
Distance<- function (par){
ysim = ode (y = y,
times = times,
func = croissance,
parms = c("Mc"= par[1],"a"= par[2],"b"=par[3]))
ssq= (yobs - ysim)^2
ssq=sum((ssq), na.rm=T) #j'ai plusieurs NA dans mon jeu de donnée d'où le na.rm = T
return(ssq) }
times = seq(from=0,to=84,by=1) #pas de temps pour l'equa diff : 1 jour
y= 2.42 #y initial pour début de calcul (valeur de la première mesure)
#optimisation
restemp <- optim(val.ini, Distance,control = list(trace=2) , method = "Nelder-Mead");restemp
je me retrouve face à 2 problèmes :
1/
D'abord l'optimisation ne fonctionne pas dans la mesure où la valeur de Mc (valeur transitoire entre les deux taux de croissance) est toujours la plus basse possible et les données collent très mal avec mon optimisation... (j'avais réalisé une optimisation manuelle des moindre carré avec des boucle for() et cela fonctionnait nettement mieux, mais les valeurs sont trop approximative, je dois donc réussir à utiliser optim()
Celà vient peut-etre de la condition if else? mais ma fonction de base (croissance) fonctionnait lorsque l'optimisation était manuelle.
2/
Ensuite si je change une valeur initiale de paramètre de Mc à 2, 3 ou + la fonction ode() bloque... la seule solution est d'augmenter le nombre de maxstep à 500 000 ou 1 000 000 ce qui fais un temps de calcul bien trop long, et je ne comprends pas d'où vient le problème
#vous pouvez essayer en mettant par exemple : val.ini=c(3, 0.01,0.01). vous pourrez voir le message d'erreur
Si vous avez des idées/suggestions qui puisse m'aider je suis preneur et vous remercie par avance.
Bien à vous,
Sylvain Bart (doctorant)