p-value

[nathalie DURAND]

Bonjour, j'ai crèé un estimateur "alpha" que je veux montrer qu'il converge vers la distribution gamma de paramètres m et 1.
pour cela j'utilise le test de kolmogorov-smirnov et comme je suis debutante pour ce logiciel R j'ai pas mal de pblm ce qui m'a poussée à venir demander de l'aide auprès de vous.
j'ai fait le code suivant qui a l'aire d'être dans le bon sens(????????????)

Code : Tout sélectionner

k<-5
n<-50
r<-40
m<-1
l<-cumsum(rexp(n)/k)
L<-exp(l)
lambda<-log((L[n]-L[n-k])/(L[n-k]-L[n-2*k]))
alpha<-k*((log(L[r+m])-log(L[r]))/lambda)
test<-ks.test(alpha,"pgamma",m,1)


mais mon problème c'est que je dois varier les k de 1 à 10 et rassembler les p-values pour tous les k affin d'en faire une graphique dont les axes sont les k et l'intervalle [0,1]
si quelqu'un peut m'aider

[Logez Maxime]

Bonjour,

je te donne juste le code pour faire ce que tu souhaites faire je ne commenterai pas la partie stat de ton analyse :

Code : Tout sélectionner

alpha <- function(k){
  n<-50
  r<-40
  m<-1
  l<-cumsum(rexp(n)/k)
  L<-exp(l)
  lambda<-log((L[n]-L[n-k])/(L[n-k]-L[n-2*k]))
  alpha<-k*((log(L[r+m])-log(L[r]))/lambda)
  test<-ks.test(alpha,"pgamma",m,1)$p.value
  test
  }

plot(1:10,sapply(1:10,alpha),ylim=c(0,1))


Maxime

[nathalie DURAND]

Merci Maxim,
mais j'ai pas bien saisi l'astuce

Code : Tout sélectionner

$p.value

dans la ligne

Code : Tout sélectionner

test<-ks.test(alpha,"pgamma",m,1)$p.value

qui fait extraire seulement les p.valeurs!
En plus si je fais sotir la fonction sapply pour voir les resultats numeriquement je remarque qu'il ne représente pas la courbe affichée et lorsque je les represente graphiquement je tombe sur deux courbe differentes alors que "normalement" doivent être la même

Code : Tout sélectionner

alpha<-function(k){
n<-50
r<-5
m<-1
l<-cumsum(rexp(n)/k)
L<-exp(l)
lambda<-log((L[n]-L[n-k])/(L[n-k]-L[n-2*k]))
alpha<-k*((log(L[r+m])-log(L[r]))/lambda)
test<-ks.test(alpha,"pgamma",m,1)$p.value
test
}
v<-sapply(1:5,alpha)
plot(1:5,sapply(1:5,alpha),type="l",ylim=c(0,1),col="blue")
lines(1:5,v,lty=1,col="red")


si vous pouvez m'expliquer ça

[Logez Maxime]

Excuse moi si je n'ai pas compris mais le but était bien de représenter les p-values obtenues pour un k variant de 1:10 en fonction de k ?

mais mon problème c'est que je dois varier les k de 1 à 10 et rassembler les p-values pour tous les k affin d'en faire une graphique dont les axes sont les k et l'intervalle [0,1]
si quelqu'un peut m'aider

Le dollars permet de récupérer que la partie de l'objet qui t'intéresse. Quand tu fais un test il te génère plein 'autres éléments que la p-value :

Code : Tout sélectionner

str(ks.test(0.5,"pgamma",1,1))
List of 5
 $ statistic  : Named num 0.607
  ..- attr(*, "names")= chr "D"
 $ p.value    : num 0.787
 $ alternative: chr "two-sided"
 $ method     : chr "One-sample Kolmogorov-Smirnov test"
 $ data.name  : chr "0.5"
 - attr(*, "class")= chr "htest"

Mais qui ne semblait pas intéressant vu ce que tu cherchais à faire. Je te laisse regarder l'aide de la fonction "$" : ?'$'.

La fonction alpha te permet juste de récupérer la p-value associé a k. Il est normal que les deux courbes ne soient pas les mêmes car dans la fonction alpha tu fias appel a un rexp. Or d'une fois sur l'autre ça ne sera pas la même séquence qui sera générée par cette fonction, donc les résultats seront différents à chaque fois que tu feras une simulation.

Si tu fais rexp(50) et si tu relances la même fonction 1000 fois tu auras autant de tirages différents.

Maxime

[nathalie DURAND]

merci ,
mais là je relance une seule fois et je trouve deux courbes differentes

[Logez Maxime]

Re,

Bien lire les précédentes réponses :

Il est normal que les deux courbes ne soient pas les mêmes car dans la fonction alpha tu fias appel a un rexp. Or d'une fois sur l'autre ça ne sera pas la même séquence qui sera générée par cette fonction, donc les résultats seront différents à chaque fois que tu feras une simulation.

rexp génère une séquence aléatoire de nombre a partir d'une loi exponentielle. Donc d'une fois sur l'autre la séquence générée ne sera pas la même et heureusement d'ailleurs car sinon pour les gens qui font des simulations ils auraient à chaque fois la même séquence de générée !

La seule façon de lui faire générer la même séquence aléatoire est de spécifier un set.seed avant le rexp. L'autre manière est de générer la séquence aléatoire avant les autres calculs, de la stocker dans un objet et de toujours faire les calculs sur cette même séquence :

Code : Tout sélectionner

alpha<-function(k,seqexp,n){
r<-5
m<-1
l<-cumsum(seqexp/k)
L<-exp(l)
lambda<-log((L[n]-L[n-k])/(L[n-k]-L[n-2*k]))
alpha<-k*((log(L[r+m])-log(L[r]))/lambda)
test<-ks.test(alpha,"pgamma",m,1)$p.value
test
}

n <- 50
sequence <- rexp(n)
sapply(1:5,alpha,seqexp=sequence,n=n)

Sachant que si tu refais sequence <- rexp(n) tu auras d'autres résultats !

Maxime

[nathalie DURAND]

c'est clair maintenant.
une dernière question ,si je veux faire cette simulation N fois pour chaque k et après je prends la moyenne pour ensuite trouver des p.values plus credibles .Que dois-je ajouter à ce code?

[Logez Maxime]

re,

Code : Tout sélectionner

alpha.sim <- function(k,N){
  # N le nombre simulations
  auxi <- function(k,n,r,m){
    l<-cumsum(rexp(n)/k)
    L<-exp(l)
    lambda<-log((L[n]-L[n-k])/(L[n-k]-L[n-2*k]))
    alpha<-k*((log(L[r+m])-log(L[r]))/lambda)
    test<-ks.test(alpha,"pgamma",m,1)$p.value
    test
    }
  n<-50
  r<-40
  m<-1
  pval <- replicate(N,auxi(k,n,r,m))
  pval <- mean(pval)
  pval
  }

plot(1:10,sapply(1:10,alpha.sim,N=1000),ylim=c(0,1))

Maxime

[nathalie DURAND]

Problème résolu.
merci

[Serge Rapenne]

Je rebondis sur ce post pour poser une question plus générale. Je débute en R mais la question d'origine m'avais semblé à ma porté. J'ai essayé de résoudre le problème et écrit du code qui fonctionne mais après avoir vu la réponse de Maxime je me rends compte que mes réflexes de programmation ne sont pas adaptés à R. J'avais en effet utilisé une boucle "for" alors que le "sapply" de Maxime me semble beaucoup plus élégant. Pouvez vous me conseiller de la doc sur l'utilisation des fonctions apply, sapply et consort ? J'ai lu l'aide en ligne mais elle ne donne aucune aide sur quand les utiliser.

Merci d'avance

Serge

[nathalie DURAND]

Salut,
voici un document qui peut t'aider surtout à partir de la page 59
http://cran.r-project.org/doc/contrib/Goulet_introduction_programmation_S.pdf
Nathalie

[nathalie DURAND]

C'est juste que j'ai oublié de mensionner que le document à l'origine c'est pour S-plus mais comme R n'est que la version gratuite de ce logiciel à mon sens avec des petites differences. Donc il peut t'aider comme il m'a aidée moi aussi sans oublier l'aide considerable et rapide des forumeurs surtout Maxime et d'autres.


Nathalie

[Pierre Bady]

bonjour

c'est pour S-plus mais comme R n'est que la version gratuite de ce logiciel à mon sens avec des petites differences.

o_O"


R n'est pas qu'une version gratuite de S-plus :')

pour information, voici "Une brève histoire de S et de R" qui nous est racontée par Jean Lobry:

http://pbil.univ-lyon1.fr/R/pdf/lang01.pdf


@+

pierre

[nathalie DURAND]

désolée de dire ça!!!c'etait pour fixer les idées.
n'empeche que c'est la realité même si R a dépassé aujourd'hui S-plus mais ses concepteurs l'ont conçu sur une base de S-plus.
lis ton document: pages 14,15,...



Nathalie

[Renaud Lancelot]

Non, non, non: la base est S, pas S+. Toutes les fonctions de base de R ont été ré-écrites "from scratch", à partir de la description de leurs syntaxe, effets et valeurs retournées telle qu'exposée dans les livres sur le langage S (*). Si vous lisez ces livres, vous trouverez d'ailleurs des fonctions qui n'ont jamais été développées ni dans R, ni dans S+. L'objectif n'est pas de "cloner" S+, mais d'assurer la portabilité entre systèmes et logiciels, ce qui est très louable. Un des exemples les plus frappants est le package lattice correspondant à la bibliothèque Trellis de S+. Le package lattice a été développé par D Sarkar à l'aide des fonctions du package grid de P. Murrell, qui n'a pas d'équivalent en S+ (grid, pas P. Murrell ;-)). Les deux package et bibliothèque ont été développés à partir de la description des fonctions (qui n'étaient alors que des maquettes) trouvées dans le livre de Bill Cleveland:

Cleveland, W.S. 1993. Visualizing data, Hobart Press.

Renaud

(*) entre autres:

(1) Richard A. Becker, John M. Chambers and Allan R. Wilks (1988), The New S Language. Chapman & Hall, New York.

(2) Chambers, J. 1998. Programming with data: a guide to the S language, Springer Verlag.