Ajustement de données

[julien hamonier]

Bonjour,

Je m'intéresse à la dissolution d'un film d'épaisseur L dans une solution aqueuse. Je recherche le coefficient de diffusion D
J'ai réalisé ce code naïf pour ajuster le paramètre D par rapport à mes données sans succès.
J'ai voulu minimiser l'erreur quadratique.

Code : Tout sélectionner

# Mes données
t <- c(0.0000, 0.0833, 0.1660, 0.2500,
       0.3330, 0.4166, 0.5000, 0.6660,
       0.8333, 1.0000, 1.5000, 2.0000,
       2.5000, 3.0000, 3.5000, 4.0000,
       4.5000, 5.0000, 5.5000, 6.0000,
       6.5000, 7.0000, 7.5000, 8.0000,
       24.0000)
y <- c(0.0000000, 0.1891365, 0.2807108,
       0.3522929, 0.3810654, 0.4089780,
       0.4273398, 0.4624385, 0.5069652,
       0.5271696, 0.5899072, 0.6566372,
       0.6942848, 0.7397635, 0.7518605,
       0.8004716, 0.8110638, 0.8310532,
       0.8723246, 0.8702028, 0.8864763,
       0.9042854, 0.9139869, 0.9285407,
       1.0000000)

data <- data.frame(t = t, y = y)

# Ma fonction d'ajustement
f <- function(t, D) {
  nbtermessomme <- 1000
  n <- 0:nbtermessomme
  L <- 0.0126
  terms <- exp(-D * (2 * n + 1)^2 * pi^2 * t / L^2) / (2 * n + 1)^2
  y <- 1 - 8 / pi^2 * sum(terms)
  return(y)
}

Dseq = seq(0,1e-4,1e-10)
res=0*Dseq
for (index in 1:length(Dseq))
   D_value <- Dseq[index]
   result <- sapply(t, function(t) f(t, D_value))
   datatemp=cbind(data,result)
   res[index]=mean((datatemp$y-datatemp$result)^2)

Dseq[which.min(res)]

Je trouve D=0 or D désigne un paramètre de diffusion et je suis surpris de trouver cette valeur. De plus, y a-t-il une méthode moins naïve pour estimer D ?
Quelqu'un pourrait-il m'aider, s'il vous plaît ?

[Logez Maxime]

Bonjour,

Il semble qu'il te manque une accolade pour ta boucle :

Code : Tout sélectionner

for (index in 1:length(Dseq)) {
   D_value <- Dseq[index]
   result <- sapply(t, function(t) f(t, D_value))
   datatemp=cbind(data,result)
   res[index]=mean((datatemp$y-datatemp$result)^2)
   }


Si tu ne la mets pas seule la première ligne sera exécutée dans la boucle :

Code : Tout sélectionner

for (index in 1:length(Dseq))
   D_value <- Dseq[index]


D'où le fait que tu ais que des 0.
Sinon tu peux utiliser des fonctions de minimisation pour obtenir ton paramètre (optim, nlm, nlminb) :

Code : Tout sélectionner

fun <- function(x) { # x = D
  result <- sapply(t, function(t) f(t, x))
  mean((data$y-result)^2)
  }
 
 nlminb(1e-10, fun, lower = 1e-12, upper = 1e-4)
$par
[1] 6.850831e-06

$objective
[1] 0.004463327


Soit un D à 6.85.10-6 ce qui est celui que tu obtiens avec ta boucle.

Cordialement,
Maxime

[julien hamonier]

Bonjour

Je vous remercie pour ce retour. J'ai manqué de vigilance sur les accolades.

J'ai également voulu comparer avec une approche de type maximum de vraisemblance

Code : Tout sélectionner

y <- c(0.0000000, 0.1891365, 0.2807108,
       0.3522929, 0.3810654, 0.4089780,
       0.4273398, 0.4624385, 0.5069652,
       0.5271696, 0.5899072, 0.6566372,
       0.6942848, 0.7397635, 0.7518605,
       0.8004716, 0.8110638, 0.8310532,
       0.8723246, 0.8702028, 0.8864763,
       0.9042854, 0.9139869, 0.9285407,
       1.0000000)
x <-diff(y)

densitef <- function(t, D) {
  nbtermessomme <- 1000
  n <- 0:nbtermessomme
  L <- 0.0126
  terms <- exp(-D * (2 * n + 1)^2 * pi^2 * t / L^2)
  y <- 8 * D / L^2 * sum(terms)
  return(y)
}


MlogL <- function(D){
  result <- sapply(x, function(t) densitef(t, D))
  logL=-sum(log(result))
  return(logL)
}
 
opt<-optimize(MlogL,c(0,1),tol=1e-20)

Je ne sais pas si vous pouvez également m'aider.
J'ai un message d'erreur : Avis : NA / Inf remplacé par la valeur maximale positive

Cordialement,
Julien

[Logez Maxime]

Re,

je pense qu'il y a un problème d'indexation dans la fonction MlogL et la valeur maximale dans optimize :

Code : Tout sélectionner

MlogL <- function(x){
  result <- sapply(t, function(u) densitef(u, x))
  logL=-sum(log(result))
  return(logL)
}
opt<-optimize(MlogL,c(0,1e-4),tol=1e-20)


Cordialement,
Maxime

[julien hamonier]

Bonjour

Je vous remercie pour votre aide.
J'ai commencé aussi à regarder du côté de la fonction nls afin de comparer mais je tombe sur un message d'erreur.
La matrice de gradient est singulière pour la valeur d'initialisation du paramètre D. Savez-vous comment choisir judicieusement une valeur initiale du paramètre D ? autrement qu'en tâtonnant . J'ai même repris la valeur de D précédemment obtenue (3.577138e-06) sans succès

Code : Tout sélectionner

## NLS
# Mes données
t <- c(0.0000, 0.0833, 0.1660, 0.2500,
       0.3330, 0.4166, 0.5000, 0.6660,
       0.8333, 1.0000, 1.5000, 2.0000,
       2.5000, 3.0000, 3.5000, 4.0000,
       4.5000, 5.0000, 5.5000, 6.0000,
       6.5000, 7.0000, 7.5000, 8.0000,
       24.0000)
y <- c(0.0000000, 0.1891365, 0.2807108,
       0.3522929, 0.3810654, 0.4089780,
       0.4273398, 0.4624385, 0.5069652,
       0.5271696, 0.5899072, 0.6566372,
       0.6942848, 0.7397635, 0.7518605,
       0.8004716, 0.8110638, 0.8310532,
       0.8723246, 0.8702028, 0.8864763,
       0.9042854, 0.9139869, 0.9285407,
       1.0000000)

mydata <- data.frame(time = t, release = y)

# Ma fonction d'ajustement
myf <- function(t, D) {
  nbtermessomme <- 1000
  n <- 0:nbtermessomme
  L <- 0.0126
  terms <- exp(-D * (2 * n + 1)^2 * pi^2 * t / L^2) / (2 * n + 1)^2
  y <- 1 - 8 / pi^2 * sum(terms)
  return(y)
}


mod1 <- nls(release ~ myf(time,D),data=mydata,start = list( D=1 ), trace=T )


Cordialement,
Julien

[Logez Maxime]

Bonjour,

Le code renvoie un warning qui dit que les longueurs d'objets ne sont pas les mêmes.
Par rapport aux codes précédents, il semble manquer un sapply dans myf.
Pour ce qui est des valeurs initiales, non je ne connais pas de technique idéale, a part essayer dans la mesure du possible que mettre une valeur dans un ordre de grandeur correcte. Ici que D soit égal à 0.1, 0.5 ou 1 ça ne change rien pour les y, probablement parce que D est bien trop élevé. Alors que lorsque D est beaucoup plus petit l'algorithme converge.
J'avoue que je ne connais pas les ordres de grandeur de D, mais avec un D de 1e-3 on y arrive, voir même avec D = 0.

Peut-être regarder du côté du package nls2 et de la fonction éponyme qui permet de tester plusieurs valeurs de paramètres de départ. Ca n'empêche qu'il faut tester des gammes de valeur appropriées.

Code : Tout sélectionner

# Ma fonction d'ajustement
myf <- function(t, D) {
  nbtermessomme <- 1000
  n <- 0:nbtermessomme
  L <- 0.0126
  y <- sapply(t, function(x) {
     terms <- exp(-D * (2 * n + 1)^2 * pi^2 * x / L^2) / (2 * n + 1)^2
     res <- 1 - 8 / pi^2 * sum(terms)
        res
       })
  return(y)
}

mod1 <- nls(release ~ myf(time,D),data=mydata,start = list( D=0), trace=T )


Cordialement,
Maxime