Tables de Mann-Whitney et Wilcoxon

[joaopa]

Bonjour, je n'ai pas trouvé les tables de Mann-Withney et de Wilcoxon sur Xcas. Je pense qu'elles devraient être impléméntées, au moins celles aux risque 5% et 1%

[parisse]

Pourquoi pas, mais n'etant pas du tout specialiste en stats, ou trouve-t-on une source expliquant comment on programme ca, si possible de maniere comprehensible pour des non specialistes en stats?

[joaopa]

Dans ce lien, il existe des méthodes (relativement facile) je crois pour implémenter la table de Mann-Whitney.
http://www.win.tue.nl/~adibucch/Nonpara ... rtmund.pdf


Je cherche pour la table de Wilcoxon

[joaopa]

Voici un programme Pari/GP qui implémente la table de Mann-Whitney avec les séries génératrices (la méthode la plus rapide à ce qu'il sembla actuellement)
Il est facilement adaptable pour en déduire la loi de distribution. Il peut être légèrement optimisé, je crois (faire une table pour les puissance par exemple au lieu de les recalculer, même si on utilise l'exponentiation rapide)
Désolé, je ne connais ni le C++ ni le python, donc je ne peux pas coder ce programme pour qu'il soit directement incorporer dans Giac

Code : Tout sélectionner

#include <stdio.h>
#include <pari/pari.h>

#define table_max 20

GEN P[table_max + 1];
GEN Q[table_max + 1][table_max + 1];
GEN R[2 * table_max + 2];
GEN S[table_max + 1][table_max + 1];
GEN U[(table_max + 1) * (table_max + 1)];

void table_P(void)
{
    int i, j;
    GEN T;

    P[0] = R[0] = mkpoln(1, gen_1);
    T = mkpoln(2, gen_1, gen_0);
    P[1] = gsub(R[0], T);
    
    for(i = 1; i <= table_max; i++)
    {
        R[i] = gsub(R[0], gpowgs(T, i));
        R[i + table_max] = gsub(R[0], gpowgs(T, table_max + i));
        P[i] = gmul(P[i-1], R[i]);
    }
    
    for(j = 0; j <= table_max; j++)
    {
        Q[0][j] = mkpoln(1, gen_1);
        for(i = 1; i <= table_max; i++)
        {
            Q[i][j] = gmul(R[i + j], Q[i - 1][j]);
            Q[i][j] = gdivent(Q[i][j], R[i]);
            S[i][j] = gdiv(Q[i][j], binomial(gadd(stoi(i), stoi(j)), j));
        }
    }
    return;
}

GEN f(int m, int n, int u)
{
    U[u] = (u == 0) ? truecoeff(S[m][n], 0): gadd(U[u-1], truecoeff(S[m][n], u));
    return U[u];
}

int main(void)
{
    int m, n, u;
    GEN alpha;

    pari_init(1000000000, 2);
    table_P();
    alpha = gdiv(stoi(5), stoi(100));

    for(m = 2; m <= table_max; m++)
        for(n = 2; n <= table_max; n++)
        {
            u = -1;
            do
            {
                u++;
            }while( (u <= m * n) && ( gcmp(gmul(f(m, n, u), gen_2), alpha) < 1) );
            if(u > 0)
                printf("%d\t", u - 1);
            else
                printf("--\t");
            if(n == table_max)
                printf("\n");
        }
        return 0;        
}

[joaopa]

D'après http://www.math.ucsd.edu/~gptesler/283/ ... andout.pdf, voici un programme pour implémenter la table de Wilcoxon

Code : Tout sélectionner

#include <stdio.h>
#include <pari/pari.h>

#define table_max 20

GEN P[table_max + 1];
GEN R[table_max + 1];
GEN U[table_max + 1];

void table_P(void)
{
    int i, j;
    GEN T;

    P[0] = R[0] = mkpoln(1, gen_1);
    T = mkpoln(2, gen_1, gen_0);
   
    for(i = 1; i <= table_max; i++)
    {
        R[i] = gadd(R[0], gpowgs(T, i));
        P[i] = gmul(P[i-1], R[i]);
        P[i] = gdiv(P[i], gen_2);
    }
   
    return;
}

GEN f(int m, int u)
{
    U[u] = (u == 0) ? truecoeff(P[m], 0): gadd(U[u-1], truecoeff(P[m], u));
    return U[u];
}

int main(void)
{
    int m, u;
    GEN alpha;

    pari_init(1000000000, 2);
    table_P();
    alpha = gdiv(stoi(5), stoi(100));
    for(m = 2; m <= table_max; m++)
    {
        u = -1;
        do
        {
            u++;
        }while( (u <= m * (m + 1) / 2) && ( gcmp(gmul(f(m, u),gen_2), alpha) < 1));
        if(u > 0)
            printf("%d\t", u - 1);
        else
            printf("--\t");
        }
        return 0;       
}

[parisse]

Le lien est plus parlant que le code pari... J'ai pour l'instant deux fonctions, une wilcoxons prenant en argument la liste de valeurs de l'echantillon et la mediane supposee et calculant la statistique W, l'autre wilcoxonp prenant en argument la taille de l'echantillon et renvoyant la liste des coefficients du polynome de la pgf, donc la liste des P(W=n). Je peux ensuite faire une fonction wilcoxont qui prendrait en arguments l'echantillon, la mediane et le niveau de confiance pour un test 2-sided.

[joaopa]

Merci beaucoup. Avoir le test prêt à l'emploi serait un plus indéniable.

Tu prévois de faire la même chose pour Mann-Whitney?

[parisse]

Bon, j'ai fini d'ecrire le test de Wilcoxon en 1 dimension. le source est a jour. Je regarderai plus tard dans la semaine ou la semaine prochaine pour le test en 2 dimensions, qui a l'air d'etre equivalent a Mann-Whitney.

[parisse]

Question: p.13 du poly sur Mann-Whitney, c'est un polynome en fait de serie generatrice, donc ca doit etre une division polynomiale, non? Du coup je ne comprends pas bien pourquoi il y a une double boucle dans le code pari...

[joaopa]

La double itération me sert à faire à faire varier la taille des échantillons (m et n)
Je ne vois pas comment on peut se passer de cette double itération si l'on veut faire varier la paire (taille_échantillon m, taille_échantillon n)

[parisse]

Ca vient surement de ce que je n'ai pas compris votre programme en pari. Je pense que je vais partir de la relation
sum_k Pmn(k)*x^k = 1/comb(m+n,n)*prod_{i=n+1}^{m+n}(1-x^i)/prod_{j=1}^m(1-x^j)

[parisse]

J'ai mis a jour le source et les versions binaires de xcas, avec les 2 tests. La commande est wilcoxont pour les 2, le 1er argument est une liste de valeurs, le 2eme soit une mediane soit une liste de valeurs, on peut rajouter un 3eme ou 4eme argument optionnel '<' ou '>' pour un test unilateral (par defaut le test est bilateral) et un reel entre 0 et 1 (alpha=0.05 par defaut). Il y a aussi la commande wilcoxonp(n) ou wilcoxonp(m,n) qui renvoie les polynomes generateurs de la distribution et wilcoxons(liste,mediane) ou wilcoxons(liste1,liste2) qui renvoie la valeur brute du test.
Il reste a tester tout ca, avant de le mettre dans l'aide.

[joaopa]

Quand je tente de compiler les sources, j'ai l'erreur suivante

pari.cc:492:41: error: ‘mpfr_get_prec’ was not declared in this scope
prec=mpfr_get_prec(e._REALptr->inf);
^
pari.cc: In function ‘giac::gen giac::in_pari(const giac::gen&, const giac::context*)’:
pari.cc:806:30: warning: cast to pointer from integer of different size [-Wint-to-pointer-cast]
argvec[k]= (GEN) v[j].val;
^
pari.cc:828:37: warning: cast to pointer from integer of different size [-Wint-to-pointer-cast]
argvec[k]=(long int*)(pos -1);
^
pari.cc:836:34: warning: cast to pointer from integer of different size [-Wint-to-pointer-cast]
argvec[k]=(long int*) v[j].val;
^
make[2]: *** [pari.lo] Erreur 1
make[2]: quittant le répertoire « /home/david/Téléchargements/giac-1.1.2/src »
make[1]: *** [all-recursive] Erreur 1
make[1]: quittant le répertoire « /home/david/Téléchargements/giac-1.1.2 »
make: *** [all] Erreur 2


Une idée pour remédier à ce problème?

[parisse]

Oui, il faut installer mpfr ou enlever
if (e.type==_REAL)
prec=mpfr_get_prec(e._REALptr->inf);
dans static GEN real2GEN(const gen & e_,GIAC_CONTEXT)

[joaopa]

J'ai testé sur les échantillons à ma dispositions. Tous les résultats obtenus sont corrects.

Cependant d'un point de vue d'un étudiant, les informations données sont peu utiles. En effet, les tables fournies aux étudiants permettent de comparer la somme signée avec des valeurs de références. On utilise pas la p-value.