Modérateur : Groupe des modérateurs
Code : Tout sélectionner
f <- function(x) (atan(x)+pi/2)*100/pi
curve(f,-10,10)
Code : Tout sélectionner
logit<-function(a2,b1,c1,t)(f(a2)/(1+((f(a2)/3)-1)*exp((b1-t)/c1)))
Code : Tout sélectionner
nls(y ~ logit(100, b1, c1, t), tab, start = list(b1= 10, c1 = 10))
François Bonnot a écrit :Je n'utilise pas non plus cette fonction, lui préférant nlm ou optim.
Logez Maxime a écrit :J'ai essayé la solution de François sur des données simulées et j'arrive à un problème de "gradient singulier" en utilisant la fonction f ou une autre pour ramener a1 entre 0 et 100. Peut-être est-ce du à mon exemple, mais peut-être que le cas va se représenter.
François Bonnot a écrit :Sur un jeu de données simulé dont la régression donne une valeur de a1 > 100, les méthodes (1) optim avec reparamétrage de a1, (2) optim avec a1 fixé à 100, et (3) optim avec upper=100 donnent des résultats quasi identiques (je peux poster le code si ça intéresse quelqu'un)
Code : Tout sélectionner
## données simulées
df <- structure(list(t = 0:50, y = c(3L, 3L, 2L, 2L, 4L, 2L, 2L, 3L,
8L, 3L, 5L, 7L, 6L, 8L, 11L, 4L, 8L, 5L, 10L, 8L, 12L, 16L, 11L,
21L, 16L, 23L, 23L, 20L, 20L, 19L, 22L, 31L, 32L, 27L, 38L, 30L,
38L, 39L, 42L, 42L, 47L, 48L, 56L, 47L, 67L, 63L, 67L, 65L, 54L,
72L, 71L)), .Names = c("t", "y"), row.names = c(NA, -51L), class = "data.frame")
plot(df,xlim=c(0,100),ylim=c(0,110))
## fonction de reparamétrage
f <- function(x) (atan(x)+pi/2)*100/pi
## fonctions correspondant ux 3 méthodes
logit1 <- function(a2,b1,c1,t)(f(a2)/(1+((f(a2)/3)-1)*exp((b1-t)/c1))) ## fonction reparamétrée
logit2 <- function(b1,c1,t)(100/(1+((100/3)-1)*exp((b1-t)/c1))) ## fonction bornée à 100
logit3 <- function(a1,b1,c1,t)(a1/(1+((a1/3)-1)*exp((b1-t)/c1))) ## fonction d'origine
## sommes de carrés à minimiser pur les 3 fonctions logit
s1 <- function(p,t,y) sum((y-logit1(p[1],p[2],p[3],t))^2)
s2 <- function(p,t,y) sum((y-logit2(p[1],p[2],t))^2)
s3 <- function(p,t,y) sum((y-logit3(p[1],p[2],p[3],t))^2)
## paramètres initaux (différents de ceux qui ont servi à simuler)
a2 <- 0 ; a1 <- f(a2) ; b1 <- 0 ; c1 <- 5
init1 <- c(a2,b1,c1)
init2 <- c(b1,c1)
init3 <- c(a1,b1,c1)
## t100 sera utilisé pour tracer les courbes complétées
t100 <- 0:100
## régression sans contrainte (avec fonction d'origine logit3)
z0 <- optim(init3,s3,t=df$t,y=df$y)
lines(t100,logit3(z0$par[1],z0$par[2],z0$par[3],t100),col="blue",lwd=2)
abline(h=z0$par[1],col="blue",lwd=2)
## méthode 1
z1 <- optim(init1,s1,t=df$t,y=df$y)
lines(t100,logit1(z1$par[1],z1$par[2],z1$par[3],t100),col="green",lwd=2)
abline(h=f(z1$par[1]),col="green",lwd=2)
## méthode 2
z2 <- optim(init2,s2,t=df$t,y=df$y)
lines(t100,logit2(z2$par[1],z2$par[2],t100),col="red",lwd=2,lty=2)
abline(h=100,col="red",lwd=2,lty=2)
## méthode 3
z3 <- optim(init3,s3,upper=100,t=df$t,y=df$y,method="L-BFGS-B") ## warning sans importance
lines(t100,logit3(z3$par[1],z3$par[2],z3$par[3],t100),col="black",lwd=3,lty=3)
abline(h=z3$par[1],col="black",lwd=3,lty=3)
## régression logistique avec glm
model <- glm(y/100~t,data=df,family=binomial)
yy1 <- model$coef[1]+model$coef[2]*t100
yy <- exp(yy1)/(1+exp(yy1))*100
## (à la fin le graphe n'est plus très lisible)
points(t100,yy,pch=2)
## méhode de simulation des données
a1 <- f(8) ; b1 <- 5 ; c1 <- 10
t <- 0:50
prob <- logit3(a1,b1,c1,t)/100
df <- data.frame(t=t,y=rbinom(length(t),100,prob))
Retourner vers « Questions en cours »
Utilisateurs parcourant ce forum : Aucun utilisateur enregistré et 1 invité