Programmation sur R erreur de valeur manquante

[invite285c27fb]

Bonjour, j'ai un dm à faire sur la logiciel R. Le but de se projet est de faire des analyse et de créer un package sur R.
J'ai un problème récurrent qui intervient sur une de mes fonctions j'ai essayé plein de procédé pour le résoudre mais je n'ai pas réussi.
Le programme consiste à estimer un paramètre grâce a la méthode de newton :

Voici mon code :

e

stimationUniv3 <- function(beta, D, eps){
  X <- D[,1]
  Y <- D[,2]
  beta_0 <- beta
  beta_1 <- beta_0 - (der1(D, beta_0)/der2(D, beta_0))
  while(abs(beta_0 - beta_1) > eps){
    a <- beta_1
    if ((der2(D, beta_1) != 0)){
      beta_1 <- beta_1 - (der1(D, beta_1)/der2(D, beta_1))
      beta_0 <- a
    }
    else {
      break
    }
  }
  return(beta_1)
}


avec der1 :
der1 <- function(D, beta){
  X <- D[,1]
  Y <- D[,2]
  d1 <- 0
  for (i in 1:length(X)){
    d1 <- d1+X[i]*(Y[i]-(exp(X[i]*beta)/(1+exp(-X[i]*beta))))
  }
  return(d1)
}


der2 :

der2 <- function(D,beta){
  X <- D[,1]
  Y <- D[,2]
  d2 <- 0
  for (i in 1:length(X)){
    d2 <- d2-((X[i]**2)*(exp(X[i]*beta)/(1+exp(X[i]*beta))**2))
  }
  return(d2)
}

D est l'échantillon fournit grâce a la fonction simulation :
simulationUniv <- function(n, beta){
  X<- rnorm(n,0,1)
  Y<-numeric(n)
  Y <-rbinom(n,1,prob=exp(t(X)*beta)/(1+(exp(t(X)*beta))))
  ech <- cbind(X, Y)
  return (ech)
}

Lorsque j'execute le programme estimationUniv l'erreur est la suivante:
Error in if ((der2(D, beta_1) != 0)) { :
missing value where TRUE/FALSE needed

cependant lorsque je teste der2 j'obtiens une valeur numérique donc je ne comprends pas comment résoudre le problème des valeurs manquantes

Merci pour votre aide
-----

[invite9dc7b526]

le message d'erreur signifie que la fonction der2 n'a pas renvoyé un nombre, mais quelque-chose qui n'est pas comparable à zéro, probablement l'une des valeurs NA (manquant) ou NaN (not a number).

[invite285c27fb]

Le problème c'est que quand je teste la fonction der2 a part j'obtiens bien un nombre.

[invite9dc7b526]

tu ne dois pas la tester avec les même arguments que quand elle est appelée par la fonction stimationUniv3.

[A voir en vidéo sur Futura]

[invite9dc7b526]

regarde:

> der2(matrix(rep(1000,4),2),1)
[1] NaN

[invite285c27fb]

D'accord dans ce cas je vais essayer de tester avec les mêmes paramètres merci

[invite285c27fb]

Ah oui d'accord mais à quoi est du ce résultat car avec d3s test sur un échantillon avec simulation celle ci me renvoyait bien un nombre

[invite6486d7bd]

Sans rien connaitre à la méthode de Newton, un truc (éventuellement) suspect (sait-on jamais) :

  d1 <- d1+X[i]*(Y[i]-(exp(X[i]*beta)/(1+exp(-X[i]*beta))))
  d2 <- d2-((X[i]**2)*(exp(X[i]*beta)/(1+exp(X[i]*beta))**2))

Vous êtes sûr que c'est un * dans la formule avec d2 ?
Ca pourrait expliquer un dépassement de capacité pour d2.

[invite285c27fb]

Oui c'est bien une multiplication, der2 et la dérivé seconde c'est enfaite la dérivé de der1 par rapport à x

[invite9dc7b526]

au fait tu n'as pas besoin de faire une boucle ici:

  d2 <- 0
  for (i in 1:length(X)){
    d2 <- d2-((X[i]**2)*(exp(X[i]*beta)/(1+exp(X[i]*beta))**2))
  }

tu peux faire juste

d2<- -sum((X^2*exp(X*beta)/(1+exp(X*beta))^2)

la fonction sum fait les itérations pour toi.

mais ça ne résoudra pas ton problème, qui est que tu calcules une exponentielle avec un argument trop grand (on est limité à exp(709) environ)

[invite6c250b59]

Les débordements viennent souvent d'une division par quelque chose trop proche de zéro. Ajouter un epsilon aux dénominateurs (genre 10^-8) pourrait résoudre le problème.

[invite285c27fb]

Bonjour, j'ai essayé en rajoutant 10^-8 au dénominateur mais cela n'a rien changé du coup j'ai remplacé le condition if par abs(der2(D, beta_1)) >0 ce qui revient au même que der2(D, beta_1)!=0 il me semble mais mon résultat est faussé, le programme n'execute pas la boucle while pourtant, quand je teste avec certaines valeurs la condition est bien vérifiée

[invite9dc7b526]

regarde aussi la fonction R all.equal()