Besoin d'aide pour un script dans nls

[Coralie Sellos]

Bonjour,

Voila je suis en stage et je dois trouver un moyen d'inclure un script dans la fonction nls (en particulier lui donner un paramètre sous forme de deux vecteurs)
Je ne vois pas du tout comment réaliser ce travail
Si l'un d'entre vous avait un exemple simple afin de me montrer comment réaliser ca, ca serait top

Merci

[marilyne berthaud]

hum....ca m'a l'air fort sympathique....
enfin avec mes stats spatiales j'ai rien à dire ..;

[Pierre Bady]

bonjour,

Voila je suis en stage et je dois trouver un moyen d'inclure un script dans la fonction nls (en particulier lui donner un paramètre sous forme de deux vecteurs)

je n'ai pas compris ce que tu voulais faire :(
la fonction 'nls' utilise un objet de type formula (cf la documentation des fonctions nls et formula).

petit exemple avec données simulées:

Code : Tout sélectionner

rndvalues <- 21
set.seed(rndvalues)
x1 <- rnorm(20,5,6)
set.seed(rndvalues^2)
x2 <- rnorm(20,9,4)
set.seed(rndvalues^3)
y <- 3+7*x1/exp(2*x2)+rnorm(20)
nls1 <- nls(y~ a+b*x1/exp(c*x2),start=list(a=1,b=6,c=3))
nls1


voir aussi:
=> http://cran.r-project.org/doc/contrib/F ... ession.pdf
=> http://pbil.univ-lyon1.fr/R/html/tdr4
=> http://www.math.unb.ca/~rolf/R/library/ ... l/nls.html
=> etc ...


HTH

pierre

[Coralie Sellos]

oui je sais, moi aussi je connais bien cette fonction
mais en fait en gros il faut que je fasse une régression avec nls mais nls doit minimiser le log de la vraisemblance pour estimer les parametres.
Je vois pas trop comment je pourrais faire
ou alors il faut que je crée une fonction qui va comprendre le nls apres ?

C'est le boulot que l'on ma donné ... mais je vois pas trop comment réaliser ca

[Logez Maxime]

Bonjour,

en espérant répondre a ta question :

Code : Tout sélectionner

x <- rnorm(100)
y <- 3*x+2+rnorm(100)
X <- data.frame(resp=y,exp=x)
file <- "nls.txt"

ted <- function(resp,exp,c,d){
  source(file,local=TRUE) # chargement du script avec ici a <- 10 et b <- 5
  pred <- c+d*exp+a-b
  (resp-pred)
  }


nls(~ted(resp,exp,c,d),data=X,start= list(c=0,d=0))
Nonlinear regression model
  model:  0 ~ ted(resp, exp, c, d)
   data:  X
     c      d
-2.999  3.017
 residual sum-of-squares: 93.08

lm(I(resp-5)~exp,data=X)
Call:
lm(formula = I(resp - 5) ~ exp, data = X)

Coefficients:
(Intercept)          exp 
     -2.999        3.017


Maxime

[Pierre Bady]

rebonjour

oui je sais, moi aussi je connais bien cette fonction
mais en fait en gros il faut que je fasse une régression avec nls mais nls doit minimiser le log de la vraisemblance pour estimer les parametres.

Dans ce cas, il me semble que les fonctions nlm, nlminb et/ou optim sont nettement plus appropriées pour ce type de problème.


HTH

pierre

[Coralie Sellos]

Oui c'est ce que je pensais, mais apparemment les hypothèses qu'il y a derrière ces fonctions sont trop restrictives pour traiter les données que l'on a.

Merci de votre aide en tout cas

[Pierre Bady]

rebonjour,

Juste par curiosité, en quoi les fonctions nlm & Co sont plus restrictives que la function nls ?


@+

pierre

[Logez Maxime]

re,

je suis aussi curieux que Pierre, et je ne vois pas très bien comment tu pourrais minimiser la vraisemblance avec nls a part s'il s'agit d'une régression linéaire classique.

Maxime

[Coralie Sellos]

Je n'en ai pas la moindre idée a vrai dire, je suis aussi sceptique que vous sur ce sujet mais bon apparemment c'est une 'ruse' qu'il a déja utlisée avec un autre logiciel... et je suis censée le faire sur R maintenant !
Désolée de ne pas pouvoir vous éclairer, j'aurais bien aimé comprendre également

Mais pour en revenir a ma question, est-il possible (par un moyen ou un autre) d'intégrer un petit script a la fonction nls ?

[Logez Maxime]

Pour en revenir à la question oui c'est possible regarde mon exemple ...

Maxime

[Pierre Bady]

rebonjour,

Mais pour en revenir a ma question, est-il possible (par un moyen ou un autre) d'intégrer un petit script a la fonction nls ?

ben désolé, mais je ne comprends toujours pas la question ...
c'est pas clair, c'est un script pour faire quoi ?

je pense que l'on a un petit problème de vocabulaire ...
pour vous, c'est quoi un script ? une function ? une formule ?

Grosso merdo, pour moi, un script est un fichier qui contient une suite de ligne de commande. exemple:

Code : Tout sélectionner

require(ade4)
?dudi.pca
data(deug)
deug.dudi <- dudi.pca(deug$tab, center = deug$cent, scale = FALSE, scan = FALSE)
par(mfrow = c(2,2))
s.class(deug.dudi$li, deug$result, cpoint = 1)

une fonction est un ensemble de commande qui réalise une tache définie (cf l'utilisation du terme 'function' dans la programmation R).
exemple:

f1 <- function(x) sum((x-mean(x))^2)/length(x)

et une formule est un objet de type 'formula' (cf la documentation):

Code : Tout sélectionner

as.formula("y~b*var1+a*exp(var2)")

Par contre, pour les nls on peut utiliser des classes spéciales 'selfStart' (elle contiennent les formules, des informations sur le calcul des valeurs initiales etc.), c'est peut-être cela qui vous intéresse.

sinon, les 'ruses' d'un soft donné ne sont pas nécessairement transposables à un autre soft ... sauf si vous avez une raison théorique ou algorithmique, mais dans ce cas-là, c'est plus une ruse.


HTH

pierre

[Coralie Sellos]

Ben c'est bien un script qu'il me faut ... Je comprend bien que vous ne compreniez pas ma question, moi même je ne comprend pas pourquoi ni comment faire ce qu'on me demande.

En gros : il faut que le nls me renvoie les paramètres d'une régression (non linéaire et avec hétéroscédasticité) mais je veux qu'il pminimise la vraisemblance( en fesant : (sqrt(Ln L)-0)² )
C'est un grand mystère pour moi mais apparemment c'est ce que je dois faire ...

Si ca n'est pas possible c'est donc normal que je n'y arrive pas :D

[Logez Maxime]

Re,

encore une fois nlm parait plus adapté (je sais je suis têtu). Voila un script pour estimer les paramètres d'un glm avec une loi de poisson avec nlm :

Code : Tout sélectionner

 #### la "-" log vraissemblance
# y la réponse
# X ton tableau des variables environnementales que tu veux utiliser
# betas qui sera estimé est le vecteur des paramètres

lnvrai3 <- function(y,X,betas){
  betas <- matrix(betas,nrow=1)
  X <- cbind(1,X)
  X <- t(X)
  if (nrow(X)!=ncol(betas))
    stop("pas le même nombre de variables que de paramètres")
  -1*(sum(-exp(betas%*%X))+sum(y*(betas%*%X))-sum(log(factorial(y))))
  }

x <- rnorm(100)
test <- rpois(100,exp(3+0.025*x))

# estimation des paramètres :
nlm1 <- nlm(lnvrai3,c(1.1,1.1),y=test,X=x,gradtol=1e-8,hessian=TRUE)

nlm1
$minimum
[1] 285.949

$estimate
[1] 2.98792327 0.01163476

$gradient
[1]  1.598049e-05 -1.273293e-05

$hessian
          [,1]      [,2]
[1,] 1987.5909  231.2939
[2,]  231.2939 2000.6188

$code
[1] 2

$iterations
[1] 9


# comparaison avec un glm

glm1 <- glm(test~x,family=poisson())

glm1

Call:  glm(formula = test ~ x, family = poisson())

Coefficients:
(Intercept)            x    #### très proche des estimations de nlm1
    2.98792      0.01164 

Degrees of Freedom: 99 Total (i.e. Null);  98 Residual
Null Deviance:      90.92
Residual Deviance: 90.65        AIC: 575.9

logLik(glm1)
'log Lik.' -285.949 (df=2) # tu retrouves le minimum de nlm1


Après a toi d'adapter ce script, tu remplaces la fonction de log vraisemblance et tu estimes tes paramètres.

Maxime

[Coralie Sellos]

Merci beaucoup je vais essayer ca ...