[Résolu] Je n'arrive pas à faire fitter mes données

Bonjour,

J'ai plusieurs jeux de données que j'aimerais faire fitter sur un modèle.
Il s'agit d'un effet sur le temps.

Mon modèle a du mal à fitter lorsqu'il y a peu de points (temps de mesure).

J'arrive à le faire fitter sur 21 points, sur 6 points également mais pas sur 4 points :

Code:


###Fit our scenarios
#21 time points
par(mfrow=c(2,1))
d=0:20
s0=52
y=c(52,62.5,65.5,66.4,65.5,64.5,63.5,62.7,61.7,60.7,59.7,58.7,57.7,56.9,55.9,54.9,54.1,53.5,52.5,52.2,52)
plot(d,y,ylab="",xlab="",ylim=c(50,70));abline(h=66.4)
a_start=1.02
Ssym_start=15
r_start=0.5
m<-nls(y~52-(a*d-Ssym*(1-exp(-r*d))),start=list(a=a_start,Ssym=Ssym_start,r=r_start))
cor(y,predict(m))
plot(d,y,ylab="",xlab="",ylim=c(50,70));abline(h=66.4)
lines(d,predict(m),col="red",lty=2,lwd=3)

a=summary(m)[11][[1]][1]
Ssym=summary(m)[11][[1]][2]
r=summary(m)[11][[1]][3]


#6 time points
par(mfrow=c(2,1))
d=0:5
s0=52
y=c(52,58.5,62.5,64.5,65.5,66.4)
plot(d,y,ylab="",xlab="",ylim=c(50,70));abline(h=66.4)
a_start=1.02
Ssym_start=15
r_start=0.5
m<-nls(y~52-(a*d-Ssym*(1-exp(-r*d))),start=list(a=a_start,Ssym=Ssym_start,r=r_start), control = list(maxiter = 500),trace = TRUE)
cor(y,predict(m))

d=0:20
a=summary(m)[11][[1]][1]
Ssym=summary(m)[11][[1]][2]
r=summary(m)[11][[1]][3]
y=s0-(a*d-Ssym*(1-exp(-r*d)))
plot(y,ylim=c(50,70));abline(h=66.4)
lines(y,col="red",lty=2,lwd=3)


#4 time points
par(mfrow=c(2,1))
d=0:3
s0=52
y=c(52,62.5,64.5,66.4)
plot(d,y,ylab="",xlab="",ylim=c(50,70));abline(h=66.4)
a_start=1.02
Ssym_start=15
r_start=0.5
m<-nls(y~52-(a*d-Ssym*(1-exp(-r*d))),start=list(a=a_start,Ssym=Ssym_start,r=r_start), control = list(maxiter = 500),trace = TRUE)
###LE MODELE NE FIT PAS
cor(y,predict(m))

d=0:20
a=summary(m)[11][[1]][1]
Ssym=summary(m)[11][[1]][2]
r=summary(m)[11][[1]][3]
y=s0-(a*d-Ssym*(1-exp(-r*d)))
plot(y,ylim=c(50,70));abline(h=66.4)
lines(y,col="red",lty=2,lwd=3)


La fonction nls() tourne en rond et finit par bloquer sur une des itérations.

Sauriez-vous si le problème vient du nombre de points que j'utilise pour fitter mes données, de la paramétrisation de départ de mes paramètres (avant de lancer la fonction nls()) ou bien d'une autre raison que je n'aurais pas identifié ?

Malheureusement je ne peux pas théoriquement augmenter mon nombre de points (4 est le maximum).

Merci d'avance

Il semblerait que je n'ai aucun résidus et que c'est cela qui fait buguer le fit.

Le modèle ne fit pas mais j'ai tout de même une liste avec des valeurs pour les paramètres donnés (celle donnée par control=list(maxiter=500) mais je ne sais pas comment récupérer ces valeur automatiquement étant donné que le modèle ne converge pas.

Solution

La fonction nls() n'aime pas les régression sans erreurs, il faut donc ajouter un aléa (j'ai fait l'erreur de vouloir fitter des données issues de la fonction du modèle... oui, je sais...).

On peut faire :

Code:

#4 time points
par(mfrow=c(2,1))
d=0:3
s0=52
y=c(52,62.5,65.5,66.4)
y=y+rnorm(length(d), sd = 0.01)
plot(d,y,ylab="",xlab="",ylim=c(50,70));abline(h=66.4)
a_start=1
Ssym_start=20
r_start=0.25
m<-nls(y~52-(a*d-Ssym*(1-exp(-r*d))),start=list(a=a_start,Ssym=Ssym_start,r=r_start), control = list(maxiter = 500),trace = TRUE)
cor(y,predict(m))

d=0:20
a=summary(m)[11][[1]][1]
Ssym=summary(m)[11][[1]][2]
r=summary(m)[11][[1]][3]
y=s0-(a*d-Ssym*(1-exp(-r*d)))
plot(y,ylim=c(50,70));abline(h=66.4)
lines(y,col="red",lty=2,lwd=3)

Ici l'aléa est y=y+rnorm(length(d),sd=0.001)

Plus de messages d'erreur de type :
"number of iterations exceeded maximum of 100"
"singular gradient"

Bonjour, pour ceux ayant toujours des problèmes à faire fitter un modèle car les données semblent trop éloignées des modèles possibles, il y a le package nls2 qui force le fit à l'aide de la méthode du maximum de vraisemblance.
On peut proposer un panel de valeurs pour chaque paramètre et le fit se fait sur le groupe de valeurs de paramètres qui minimise le Residual Sum of Squares.

Salut,

Je ne suis pas sûr de comprendre ce que tu appelles "forcer".
Disons que nls2 permet d'éviter d'avoir à se taper "à la main" le test de nombreuses valeurs de départ des paramètres quand l'algorithme a du mal à se caler au départ.

Nik

Bonjour Nik,

En faite nls() ne va pas forcément converger, parfois s'il arrive à trouver un pool de paramètres pour lesquels il trouvera un RSS nul (Residual Sum of Square), le code semble buguer.

nls2 va forcément converger si on lui précise qu'on veut forcément trouver un pool de paramètres, c'est ça que j'entends par "forcer".

En faite j'ai l'impression qu'on applique le maximum de vraisemblance.
Penses-tu que ce soit la meilleure méthode pour obtenir les paramètres estimés d'une courbe ?