[Statistique] Test du "khideux" hors des clous

Le
Francois Grieu
Bonjour,

problème (non scolaire, voir section E.6 page 30 de [trngk31e])
où l'on fait un test du "khideux" loin du domaine d'application
habituel.

On tire N€ fois un entier supposé équidistribué de 1 à K;
on compte le nombre Mj d'occurrences de chaque valeur;
on calcule

X = Somme (Mj - N/K)^2 * K/N
j=1..K

On demande la probabilité p(65) que X>65 (événement improbable
correspondant à une distribution irrégulière des Mj).

J'ai fait une simulation, qui donne
p ~= 3,7*10^-7 (seul l'ordre de grandeur est sur)

La source du problème (section E.6 page 30 de [trngk31e])
indique
p ~= 3,8*10^-7
et explique que l'on obtient ce résultat en considérant que X
suit une distribution du "khideux" avec 15 degrés de liberté,
citant [Ka].
Pourtant, avec cette hypothèse, on obtient
p ~= 3,4*10^-8
Par exemple avec Excel =LOI.KHIDEUX(65;15)


Il me semble plausible que, compte tenu de la probabilité
très faible considérée, et/ou du petit N/K, l'hypothèse
d'une distribution du "khideux" soit complètement fausse.

Mais, mis à part une simulation, comment faire ?


Note: [Ka] limite son test du "khideux" à des probabilités
beaucoup moins petites.

François Grieu



[trngk31e] W. Killmann & W. Schindler: Functionality classes
and evaluation methodology for true (physical) random number
generators
http://www.bsi.de/zertifiz/zert/interpr/trngk31e.pdf

[Ka] G.K. Kanji: 100 Statistical Tests
Vidéos High-Tech et Jeu Vidéo
Téléchargements
Vos réponses Page 1 / 2
Gagnez chaque mois un abonnement Premium avec GNT : Inscrivez-vous !
Trier par : date / pertinence
Acetonik
Le #17697521
"Francois Grieu"
Bonjour,



[...]
On tire N€ fois un entier supposé équidistribué de 1 à K;
on compte le nombre Mj d'occurrences de chaque valeur;
on calcule

X = Somme (Mj - N/K)^2 * K/N
j=1..K



[...]
J'ai fait une simulation, qui donne
p ~= 3,7*10^-7 (seul l'ordre de grandeur est sur)

La source du problème (section E.6 page 30 de [trngk31e])
indique
p ~= 3,8*10^-7
et explique que l'on obtient ce résultat en considérant que X
suit une distribution du "khideux" avec 15 degrés de liberté,
citant [Ka].



En toute rigueur X converge vers une variable "khideux" lorsque N tend
vers l'infini , mais ne suit pas exactemement une loi du "khideux" pour N
fixé.

On admet en pratique que X suit une loi du "khideux" lorsque les Mj sont
tous supérieurs à 10, parfois même lorsque les Mj sont tous supérieurs à
5 , ce qui est le cas de cet exemple.


Pourtant, avec cette hypothèse, on obtient
p ~= 3,4*10^-8
Par exemple avec Excel =LOI.KHIDEUX(65;15)



X n'étant pas exactement une variable "khideux" , il est normal que l'on
obtienne une distorsion entre une simulation ( pseudo aléatoire) et les
résultats obtenus par l'application de la loi du "khideux"


Il me semble plausible que, compte tenu de la probabilité
très faible considérée, et/ou du petit N/K, l'hypothèse
d'une distribution du "khideux" soit complètement fausse.



C'est vrai que N/K=5 est un cas limite pour considérer que X suit
approximativement une loi du "khideux" , cependant pas tres éloigné du
domaine d'application habituel.

Mais, mis à part une simulation, comment faire ?


On doit de toute façon s'en contenter car la loi de X n'est pas connue ( Les
Mj ne sont pas indépendantes leur somme est N...dommage !!)

Cordialement
Acetonik
Francois Grieu
Le #17698311
Dans l'article "Acetonik"
"Francois Grieu"
> On tire N€ fois un entier supposé équidistribué de 1 à K;
> on compte le nombre Mj d'occurrences de chaque valeur;
> on calcule
>
> X = Somme (Mj - N/K)^2 * K/N
> j=1..K



Rappel: On demande la probabilité p(65) que X>65 (événement
improbable correspondant à une distribution irrégulière des Mj).
L'hypothèse standard d'une distribution "khideux" donne un résultat
faux par un facteur plus de 10 (3,4*10^-8 au lieu de 3,8*10^-7).

En toute rigueur X converge vers une variable "khideux" lorsque N tend
vers l'infini, mais ne suit pas exactemement une loi du "khideux" pour N
fixé.

On admet en pratique que X suit une loi du "khideux" lorsque les Mj sont
tous supérieurs à 10, parfois même lorsque les Mj sont tous supérieurs à
5 , ce qui est le cas de cet exemple.

C'est vrai que N/K=5 est un cas limite pour considérer que X suit
approximativement une loi du "khideux" , cependant pas tres éloigné du
domaine d'application habituel.



Je pense que mon problèle est la combinaison de Mj "petit" et d'une
zonne très "à droite" de la distribution, à très faible probabilité.

> Mais, mis à part une simulation, comment faire ?



On doit de toute façon s'en contenter car la loi de X n'est pas connue
( Les Mj ne sont pas indépendantes leur somme est N...dommage !!)



Pour de très petits N je parviens à calculer le distribution
exacte (sans simulation). Pour N plus grand, je sèche, mais je
ne déserpère par d'avoir la distribution exacte pour N€
et les grands X, ce qui m'intéresse


François Grieu
Acetonik
Le #17699291
"Francois Grieu"

> On tire N€ fois un entier supposé équidistribué de 1 à K;
> on compte le nombre Mj d'occurrences de chaque valeur;
> on calcule
>
> X = Somme (Mj - N/K)^2 * K/N
> j=1..K



Rappel: On demande la probabilité p(65) que X>65 (événement
improbable correspondant à une distribution irrégulière des Mj).




L'hypothèse standard d'une distribution "khideux" donne un résultat
faux par un facteur plus de 10 (3,4*10^-8 au lieu de 3,8*10^-7).


Evidemment puisque , comme déjà dit, X ne suit pas une loi du "Khi deux".
(X suit approximativement une loi du "Khi deux" seulement si N est assez
grand....avec les Npj supérieurs à 10 , le plus souvent en pratique )


[...]
Pour de très petits N je parviens à calculer le distribution
exacte (sans simulation).



Quelle est la loi de X pour N=3 par exemple?
( Les Mj suivent des loi binomiales NON indépendantes)

Pour N plus grand, je sèche, mais je
ne déserpère par d'avoir la distribution exacte pour N€


Alors là , je doute fort que vous obteniez la distribution exacte de X .

et les grands X, ce qui m'intéresse


ne veut rien dire ...


Bon courage, :-)
Francois Grieu
Le #17699721
"Acetonik"
"Francois Grieu"
On tire N€ fois un entier supposé équidistribué de 1 à K;
on compte le nombre Mj d'occurrences de chaque valeur;
on calcule

X = Somme (Mj - N/K)^2 * K/N
j=1..K

Rappel: On demande la probabilité p(65) que X>65 (événement
improbable correspondant à une distribution irrégulière des Mj).



[...]
Pour de très petits N je parviens à calculer le distribution
exacte (sans simulation).



Quelle est la loi de X pour N=3 par exemple?



Prenons N=3 et gardons K. Convenons de classer les Mj
par ordre décroissant. Il n'y a que 3 cas à considérer:
M0=3 => XE probabilité 1/256
M0=2, M1=1 => X#+2/3 probabilité 45/256
M0=1, M1=1, M2=1 => X probabilité 210/256

Je sais par exemple que X>65 avec une probabilité nulle,
et que X>15 avec une probabilité 23/128.


( Les Mj suivent des loi binomiales NON indépendantes)



Pour des petits nombres seules les lois discrètes sont exactes.

Pour N plus grand, je sèche, mais je
ne déserpère par d'avoir la distribution exacte pour N€





Alors là , je doute fort que vous obteniez la distribution exacte de X.

et les grands X, ce qui m'intéresse


ne veut rien dire ...



Précisons: pour N€, je désespère (pour l'instant) de calculer
tout le tableau, mais je peux calculer le haut, par X décroissant
M0€ => X00 probabilité 1/2^316
M0y M1=1 => X68+2/5 probabilité 240/2^316
M0x M1=2 => X37+3/5 probabilité ...
M0x M1=1 M2=1 => X37+1/5 probabilité ...

Et je sais par exemple que X>1150 avec une probabilité 241/2^316
Si je pouvais descendre à X>65, j'aurais ma réponse...


Bon courage, :-)



Merci.


François Grieu
Acetonik
Le #17703521
"Francois Grieu"
"Acetonik"



Précisons: pour N€, je désespère (pour l'instant) de calculer
tout le tableau, mais je peux calculer le haut, par X décroissant
M0€ => X00 probabilité 1/2^316
M0y M1=1 => X68+2/5 probabilité 240/2^316
M0x M1=2 => X37+3/5 probabilité ...
M0x M1=1 M2=1 => X37+1/5 probabilité ...



Ok, je comprends mieux ce que vous voulez .
Si c'est la loi conjointe des Mj que vous cherchez , pour n'obtenir que le
résultat P(X>65) dans le cas particulier où N€ , alors cela pourrait être
envisageable .

Démarche succincte:( impossible à la main , peut être programmable?)

Partir de la conjointe des Mj ( j variant de 1 à 16 avec la notation
première ) qui est une loi multinomiale , voir par exemple:
http://nte-serveur.univ-lyon1.fr/immediato/Math/Enseignement/06%20Probabilit%E9s/11.%20Loi%20multinomiale/D%E9finitions_11.htm

lien en une seule ligne: http:// ........E9finitions_11.htm


On a donc :

P[M1=m1, M2=m2 , ... , M16=m16]
N! / (m1! m2! ... m16!) * (p1^m1) * (p2^m2) *...* (p16^m16)
80! / (m1! m2!... m16!) * (1/16)^80

Retenir tous les cas [M1=m1, M2=m2 , ... , M16=m16] donnant X>65
Puis les additionner..... c'est tout !!!

Cordialement
Acétonik

PS/
Appliquer la loi du "Khideux" ,dans ce cas certes très limite, donne
finalement et très rapidement un résultat qui n'est pas si aberrant que cela
en comparant avec la simulation. (ordres de grandeur 10^-8 et 10^-7)
Acetonik
Le #17706971
j'ajoute qu'il y a 110 375 347 398 090 219 solutions entieres positives

à l équation m1+m2+......+m16 = 80

résultat obtenu par la formule classique :
http://www.bibmath.net/dico/index.php3?action¯fiche&quoi=./c/combirepet.html

Prévoyez un super calculateur pour trouver tous les cas donnant X>65
....!!!

Acetonik
Francois Grieu
Le #17723751
"Acetonik" a écrit:
j'ajoute qu'il y a   110 375 347 398 090 219  solutions entieres posi tives

à  l équation   m1+m2+......+m16 = 80

résultat obtenu par la formule classique :http://www.bibmath.net/dico/i ndex.php3?action¯fiche&quoi=./c/combir...

Prévoyez un super calculateur pour trouver tous les cas  donnant X>65
....!!!



Certes, mais
- on peut se restreindre aux solutions avec m1>=m2>=m3..>=m16, ce qui
gagne 13 ordres de grandeur;
- et on n'est intéressé que par les solutions qui donnent X>65, ce qui
gagne encore 6 ordres de grandeur

François Grieu
Acetonik
Le #17728031
"Francois Grieu"
"Acetonik" a écrit:
>j'ajoute qu'il y a 110 375 347 398 090 219 solutions entieres positives


En réalité il faut envisager encore bien plus de 16-uplets pour ne retenir
que les solutions qui conviennent même si l'on peut forcer quelques boucles
à se terminer sous certaines conditions ....

- on peut se restreindre aux solutions avec m1>=m2>=m3..>=m16, ce qui
gagne 13 ordres de grandeur;


Oui , mais il faut alors tester le nombre de répétions des mj (et ces test
prennent du temps!!) ,on obtient par exemple:
en tenant compte de l'ordre:

avec ordre w1=(28,10,10,5,5,5,5,2,2,2,1,1,1,1,1,1)
p(w1)== 80! / (28!10!10!5!5!5!5!2!2!2!1!1!1!1!1.1!)*( 1/16)^80
et
sans ordre w2={28,10,10,5,5,5,5,2,2,2,1,1,1,1,1,1}
p(w2)== 16!/ (1!2!4!3!6!)*p(w1)


- et on n'est intéressé que par les solutions qui donnent X>65, ce qui
gagne encore 6 ordres de grandeur


Oui mais il faut trier les (m1,m2,....m16) qui donnent X>65 de façon
efficace...
Toujours pour K, j'obtiens P(X>65) en valeur exacte, pour les valeurs :
N= 3,4,5,6,7,8,9,10,11,12,13 ( en moins de 2 heures pour N).
Le temps de calcul est en gros multiplié par 2 à chaque valeur suivante de N
par exemple:
N p(X>65) 887189 / 4398046511104 = 4,749196933*10^-6
N p(X>65)8098927 / 35184372088832 = 3,925007576*10^-6


Il semble bien qu'il soit impossible d'envisager d'aller jusqu'à N€...
(cela me rappelle un peu le problème des grains de blé sur l'échiquier)
Je serais d'ailleurs curieux de savoir jusqu'à quelle valeur de N vous
pouvez aller en un temps raisonnable, disons en gros 2 heures , sachant que
cela dépend bien sûr des processeurs (1733Mhz pour moi)

Cordialement
Acetonik
Francois Grieu
Le #17741271
Dans l'article "Acetonik" a écrit:

"Francois Grieu"
"Acetonik" a écrit:
j'ajoute qu'il y a 110 375 347 398 090 219 solutions entieres positives




En réalité il faut envisager encore bien plus de 16-uplets pour ne retenir
que les solutions qui conviennent même si l'on peut forcer quelques boucles
à se terminer sous certaines conditions ....



Pas compris pourquoi votre 110 375 347 398 090 219 = C(95,15)
ne suffirait pas.

- on peut se restreindre aux solutions avec M1>=M2>=M3..>=M16, ce qui
gagne 13 ordres de grandeur;





Là j'ai été nettement trop optimiste; je suis parti de log10(16!) ce
qui marcherais pour N grand et la totalité d l'exploration, mais avec
N€ il y a beaucoup de Mj égaux, et particulièrement dans la partie à
explorer :-(

Oui, mais il faut alors tester le nombre de répétions des Mj (et
ces test prennent du temps!!)
on obtient par exemple: en tenant compte de l'ordre:

avec ordre w1=(28,10,10,5,5,5,5,2,2,2,1,1,1,1,1,1)
p(w1)== 80! / (28!10!10!5!5!5!5!2!2!2!1!1!1!1!1.1!)*( 1/16)^80
et
sans ordre w2={28,10,10,5,5,5,5,2,2,2,1,1,1,1,1,1}
p(w2)== 16!/ (1!2!4!3!6!)*p(w1)



Oui. Mais on peut pré-calculer le terme au dénominateur (1!2!4!3!6!)
dans l'expression ci-dessus; il n'y en a que 2^152768, que l'on peut
indexer par l'entier formé par la suite de 15 bits définie par
(M1=M2),(M2=M3),..(M15=M16) laquelle peut être construite au fur
et à mesure de l'exploration.

- et on n'est intéressé que par les solutions qui donnent X>65, ce qui
gagne encore 6 ordres de grandeur


Oui mais il faut trier les (m1,m2,....m16) qui donnent X>65 de façon
efficace...



On peut explorer en maintenant le tri. Dans ce qui suit je considère
seulement les (M1,M2..,M16) avec M1>=M2>=M3..>=M16 et M1+M2..+M15€,
que j'explore par ordre lexicograhique inverse et en calculant X
M1 M2 M3 M4 M5 M6 M7 M8 M9 M10 M11 M12 M13 M14 M15 M16 X
80 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1200
79 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1168,4
78 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1137,6
78 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1137,2
77 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1107,6
77 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1106,8
76 4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1078,4
76 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1077,2
76 2 2 0 0 0 0 0 0 0 0 0 0 0 0 0 1076,8
76 2 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1076,4
76 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 1076
75 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1050
(..)
65 2 2 1 1 1 1 1 1 1 1 1 1 1 0 0 768,8
65 2 1 1 1 1 1 1 1 1 1 1 1 1 1 0 768,4
65 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 768
64 16 0 0 0 0 0 0 0 0 0 0 0 0 0 0 790,4
64 15 1 0 0 0 0 0 0 0 0 0 0 0 0 0 784,4
64 14 2 0 0 0 0 0 0 0 0 0 0 0 0 0 779,2
64 14 1 1 0 0 0 0 0 0 0 0 0 0 0 0 778,8
(..)
25 16 10 3 3 3 3 3 2 2 2 2 2 2 1 1 130,4
25 16 10 3 3 3 3 2 2 2 2 2 2 2 2 1 130
25 16 10 3 3 3 2 2 2 2 2 2 2 2 2 2 129,6
25 16 10 4 4 4 4 4 4 4 1 0 0 0 0 0 138,8
25 16 10 4 4 4 4 4 4 3 2 0 0 0 0 0 138
25 16 10 4 4 4 4 4 4 3 1 1 0 0 0 0 137,6
25 16 10 4 4 4 4 4 4 2 2 1 0 0 0 0 137,2
25 16 10 4 4 4 4 4 3 3 3 0 0 0 0 0 137,6
25 16 10 4 4 4 4 4 3 3 2 1 0 0 0 0 136,8
25 16 10 4 4 4 4 4 3 3 1 1 1 0 0 0 136,4
(..)

L'ordre d'exploration coïncide **grossièrement** avec l'ordre
décroissant de X, ce qui est favorable. Et quand (si jamais..) on
atteind le premier terme avec X<e il semble possible d'élaguer une
partie de l'exploration pour ne retenir que les termes suceptibles
de produire X>65.

Mais je ne parviens pas à dénombrer le nombre de pas nécessaires,
même approximativement.

Toujours pour K, j'obtiens P(X>65) en valeur exacte, pour les valeurs :
N= 3,4,5,6,7,8,9,10,11,12,13 ( en moins de 2 heures pour N).
Le temps de calcul est en gros multiplié par 2 à chaque valeur suivante de N
par exemple:
N p(X>65) 887189 / 4398046511104 = 4,749196933*10^-6
N p(X>65)8098927 / 35184372088832 = 3,925007576*10^-6


Il semble bien qu'il soit impossible d'envisager d'aller jusqu'à N€...
(cela me rappelle un peu le problème des grains de blé sur l'échiquier)
Je serais d'ailleurs curieux de savoir jusqu'à quelle valeur de N vous
pouvez aller en un temps raisonnable, disons en gros 2 heures , sachant que
cela dépend bien sûr des processeurs (1733Mhz pour moi)




Je onjecture que mon exploration ci-dessus est moins vorace que O(2^K).
On en reparle dans une semaine !


François Grieu
Publicité
Poster une réponse
Anonyme