programme tirant une boule parmi d'autres, pondéré avec des 'poids'

[gloups13]

Bonsoir,
Je rencontre actuellement un problème qui, à première vue, pourrait paraître simple mais en y réfléchissant bien, n'est pas si simple que ça.
Je vous l'explique:
J'ai une matrice B de 1 ligne et 190 colonnes.
Dans chaque colonne, il y a une probabilité. Les probabilités sont rangées par ordre croissant dans B.
J'aimerai faire un programme qui me sélectionne une colonne de B en tenant compte de la probabilité de sélectionner une colonne (probabilité qui est inscrite dans B.
J'ai essayé de modéliser les cases de B par des boules et de tirer une boule parmi les 190 en tenant compte de leurs poids. ça n'a pas marché.
J'ai essayé de convertir les probabilité en angle pour utiliser une roue comme il y a dans les casinos, ça a été un échec. Je sais vraiment pas comment faire. Si vous avez une idée.
Remarque: j'ai des probabilités très petites: ça va de 10^-8 à10^-2
J'utilise Scilab pour coder.
Merci de votre aide.
CDT
-----

[Evil.Saien]

Salut,

C'est quoi ton code actuel qui ne marche pas ?

Bon, peut-être que tu veux explorer une nouvelle piste. T'as une idée de par où commencer ?

[gloups13]

Bonsoir, alors voici le programme
Ici ce programme ne marche même pas alors que la matrice qui contient les probabilité (ici la matrice P) est de taille petite.

P=[3/6, 2/6, 1/6];//Matrice qui contient la probabilité de tirer le chromosome 1, 2, 3...N
I=[];
ADN_next=zeros(N,tc);
for i=1:N-1
    for j=i+1:N
        p=P(1,i)*P(1,j)*(1/(1-P(1,i))+1/(1-P(1,j)));//Probabilité de séléctionner les chromosomes i et j.
        teta=fix(360*p);
        B=10*(i+j/10);// Noyau du programme.ici, B n'est rien d'autre que le nombre ij (c'est pas i*j)!!!!
        I=[I,B*ones(1,teta)];//I est la roue biaisée du casino. 
    end
end

et le résultat:

         column 1 to 12

   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.

         column 13 to 24

   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.

         column 25 to 36

   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.

         column 37 to 48

   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.

         column 49 to 60

   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.

         column 61 to 72

   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.

         column 73 to 84

   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.

         column 85 to 96

   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.

         column 97 to 108

   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.

         column 109 to 120

   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.

         column 121 to 132

   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.

         column 133 to 144

   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.

         column 145 to 156

   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.

         column 157 to 168

   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.

         column 169 to 180

   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.

         column 181 to 192

   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.

         column 193 to 204

   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.   12.

         column 205 to 216

   12.   12.   12.   12.   12.   13.   13.   13.   13.   13.   13.   13.

         column 217 to 228

   13.   13.   13.   13.   13.   13.   13.   13.   13.   13.   13.   13.

         column 229 to 240

   13.   13.   13.   13.   13.   13.   13.   13.   13.   13.   13.   13.

         column 241 to 252

   13.   13.   13.   13.   13.   13.   13.   13.   13.   13.   13.   13.

         column 253 to 264

   13.   13.   13.   13.   13.   13.   13.   13.   13.   13.   13.   13.

         column 265 to 276

   13.   13.   13.   13.   13.   13.   13.   13.   13.   13.   13.   13.

         column 277 to 288

   13.   13.   13.   13.   13.   13.   13.   13.   13.   13.   13.   13.

         column 289 to 300

   13.   13.   13.   13.   13.   13.   13.   13.   13.   13.   13.   13.

         column 301 to 312

   13.   13.   13.   13.   13.   23.   23.   23.   23.   23.   23.   23.

         column 313 to 324

   23.   23.   23.   23.   23.   23.   23.   23.   23.   23.   23.   23.

         column 325 to 336

   23.   23.   23.   23.   23.   23.   23.   23.   23.   23.   23.   23.

         column 337 to 348

   23.   23.   23.   23.   23.   23.   23.   23.   23.   23.   23.   23.

         column 349 to 358

   23.   23.   23.   23.   23.   23.   23.   23.   23.   23.

remarquez que I n'a même pas 360 colonnes. C'est du au f que les p sont trop petits.

C'est difficile de comprendre ce que j'ai fait. En fait, dans I il y a un nombre (par exemple 23) cela veut dire que si en générant aléatoirement un nombre entier (compris entre 1 et 358), par exemple 350, je sélectionne la 350ème colonne de I et je relève sa valeur: ici '23' alors le chromosome 2 et 3 ont étaient sélectionnés.
Ensuite, comme j'ai dis qu'il y avait des poids, j'ai représentais ces poids aves un nombre de colonne déterminé qui comportent le même nombre.

Dans ce cas j'ai que 3 chromosomes et ça ne marche pas bien.
Hors ce cas est un exemple, dans mon projet j'ai beaucoup plus de chromosomes alors cette méthode ne marchera vraiment pas.
Je sais pas trop comment faire

[Evil.Saien]

Commençons simple.

Admettons que tu as 2 boules. La probabilité de tirer la boule A est pA et celle de la boule B est pB (=1-pA).

Comment tu ferais pour faire un tirage pondéré ?

[A voir en vidéo sur Futura]

[gloups13]

Alors soit P_A=0.7 et P_B=0.3
Je crée la matrice I=[1 1 1 1 1 1 1 2 2 2 ]
Je génère un nombre entier aléatoire entre 1 et 10
Si le nombre généré aléatoirement vaut, par exemple, 4 alors je sélectionne I(1,4) qui vaut 1 ce qui veut dire que la boule A a été tirée.
Cette méthode ne marche plus quand ya beaucoup de boules avec des probabilités très petites.

[Evil.Saien]

Ok, t'as le principe.

Maintenant essaie de le généraliser avec n'importe quels pA et pB (a priori inconnus, on sait juste que pA >= 0, pB >=0, pA+pB=1)

[gloups13]

Al là je n'y arrive pas. Car si P_A=10^-8 ,comment faire un programme qui tienne compte de ce cas. Je vois pas comment faire.

[Evil.Saien]

Oui je comprends que tu n'y arrives pas, parce que tu énumères explicitement les possibilités.

Et si tu oubliais ta matrice I pour raisonner par intervalles...

[gloups13]

comme je l'ai dit, j'ai beaucoup cherché avant de poster cette question. J'ai essayé les intervalles mais si j'ai un intervalle I_A=[0,P_A] et I_B=[P_A,1[
Si je génère un nombre aléatoire entre 0 et 1 alors je ne tiens plus compte du fait que les boules sont pondérées Car la probabilité de tirer un nombre dans l'un des deux intervalle est la même car tous les deux comportent une infinité de nombre.

[Dlzlogic]

Bonsoir,
Il y a aussi un point important qui semble oublié : la somme des probabilités des 190 boules doit être égale à 1, c'est à dire 100%.
Pour ce type de problème, j'avais conseillé la création de classes.
Encore une fois, pardon si je me répète, si on s'obnubile sur la méthode de résolution avant d'avoir posé le problème, on a peu de chance d'y arriver.

[gloups13]

Bonsoir,
Pas d'inquiétude, aussi petites quelle soit, la somme de toutes les probabilités fait bien 1
Ensuite, je ne rejette aucune méthode. Vous dites avoir évoqué la création de classe. Je vais aller revoir ce que vous avez dit sur ce sujet. Mais les classe en Scilab, je ne crois pas que ça existe.
Bref, c'est pas simple cette histoire.

[Dlzlogic]

Le terme de "classe" n'a rien à voir avec les classes des logiciels modernes.
Oubliez l'informatique, écrivez l'algorithme, c'est à dire la logique de ce que vous voulez faire, il sera bien temps, ensuite, de le traduire en code, dans votre langage préféré.

[Evil.Saien]

comme je l'ai dit, j'ai beaucoup cherché avant de poster cette question. J'ai essayé les intervalles mais si j'ai un intervalle I_A=[0,P_A] et I_B=[P_A,1[
Si je génère un nombre aléatoire entre 0 et 1 alors je ne tiens plus compte du fait que les boules sont pondérées Car la probabilité de tirer un nombre dans l'un des deux intervalle est la même car tous les deux comportent une infinité de nombre.

Ce n'est pas parce qu'il existe une bijection entre [0 1e-36] et [0.5 1] qu'il y a autant de chances qu'un nombre choisi aléatoirement entre 0 et 1 ait autant de chances de se trouver dans chacun des intervales.

Petite démonstration: prenons A un réel nombre choisi aléatoirement dans [0 1] selon une distribution uniforme. T'es d'accord pour dir que p(A<=0.5) = 0.5, et p(A>=0.5) = 0.5.
Maintenant, disons que A>0.5, selon le même principe, on a p(A>=0.75|A>=0.5) = 0.5 et p(A<=0.75|A>0.5) = 0.5.

Au final, qu'elle est la probabilité que A > 0.75 ? Selon le théorème de Bayes:
p(A>=0.75|A>=0.5)*p(A>=0.5) = p(A>=0.5|A>=0.75)*p(A>=0.75)

On voit clairement que p(A>=0.5|A>=0.75) = 1, d'où:
p(A>=0.75) = p(A>=0.75|A>=0.5)*p(A>=0.5) = 0.5*0.5 = 0.25

Donc, bien que l'interval [0.75 1] ait autant d'éléments que [0 0.75], A a 25% de chance de s'y trouver. Tu noteras que la probabilité de 0.25 est égale à la longueur de [0.75 1] et que ce n'est pas le fruit du hasard.

[Dlzlogic]

@ Evil,
C'est vraiment casse-cou ce que tu expliques.
D'abord, le seul moyen de tirage direct que je connaisse avec une machine, est le résultat de la fonction rand. C'est à dire que c'est un tirage aléatoire, uniformément distribué et qui satisfait aux lois des probabilités.
D'autre part, je ne vois pas ce que vient faire Bayes dans l'histoire. Comment utiliser "sachant que", qui sait quoi ?
Il me parait indispensable que notre ami sache que si on fait un tirage entre 0 et 1, la probabilité d'avoir un résultat de l'ordre de 0.1 est exactement la même que celle d'avoir un résultat de l'ordre de 0.9, ou d'avoir un résultat de l'ordre de 0.6.

[invite73192618]

+1

Autrement dit, tout ce que tu as besoin de faire est de creer un vecteur contenant les bornes, puis de tirer au hasard entre 0 et 1, et finalement de trouver la borne correspondante. Le tout se fait en quelques lignes.

c(1)=0; for i=1:190; c(1+i)=b(i)+c(i); end; 
disp(max(find(c<rand(1)))

A noter que, dependant de ce que tu veux faire quand ca tombe sur une borne, tu pourrais choisir min et > plutot que max et < a la derniere ligne.

[Evil.Saien]

c(1)=0; for i=1:190; c(1+i)=b(i)+c(i); end; 
disp(max(find(c<rand(1)))

Ce code est fonctionel mais pas optimisé. De manière générale, si on peut éviter une boucle (dans un language vectoriel) alors il vaut mieux l'éviter.
find(c<rand(1)) peut être assez long si c est très grand.

Ce serait beaucoup mieux optimisé en écrivant ceci (sous matlab) :

c = cumsum(P);
n_rand = find(c < rand(1), 1, 'last');
disp(num2str(n_rand));

[Evil.Saien]

D'abord, le seul moyen de tirage direct que je connaisse avec une machine, est le résultat de la fonction rand

Ah bah si c'est le seul que tu connaisses, c'est sans doute le seul qui doit exister alors...

D'autre part, je ne vois pas ce que vient faire Bayes dans l'histoire. Comment utiliser "sachant que", qui sait quoi ?

Tu verras peut-être mieux en lisant ceci:
https://fr.wikipedia.org/wiki/Th%C3%...%A8me_de_Bayes

En général, plutôt que d'affirmer "ceci est alambiqué / farfelus / faux / "casse-cou"", suffit de dire "j'ai pas compris". On t'explique et c'est plus productif

[invite73192618]

+1

Juste pour préciser: +1 message #13.

Ce code est fonctionel mais pas optimisé.

Tout-à-fait, avec quelques centaines d'éléments cela ne fera pas de différence mais ton écriture est indiscutablement meilleure.

[Dlzlogic]

Bonjour Evil,
Oui, peut-être fais-tu allusion aux fonctions mises en œuvre par des ordinateurs quantiques.
Sinon, donne moi un exemple de méthode de tirage direct différent de la méthode rand, s'il te plait.

La théorème de Bayes est basé sur le fait que lors d'évènements précédents on a modifié les éléments servant au tirage suivant, d'où l'expression "sachant que". D'où ma question : "Qui sait Quoi ?". Dans le cas du sujet dont il est question, le chromosome peut être tiré plusieurs fois de suite, mais chaque fois sa note aura changé.
En d'autres terme, on tire deux chromosomes, chacun ayant sa probabilité d'être tiré, calculé avec sa note, cela produit deux autres chromosomes avec chacun leur note, donc leur probabilité, et on recommence.

[Evil.Saien]

Tout-à-fait, avec quelques centaines d'éléments cela ne fera pas de différence mais ton écriture est indiscutablement meilleure.

T'as raison, ça changera même pas d'un iota la vitesse du code

J'avoue qu'optimiser les codes matlab (surtout enlever les boucles) fait partie de mes petits plaisirs

[Evil.Saien]

Bonjour Evil,
Oui, peut-être fais-tu allusion aux fonctions mises en œuvre par des ordinateurs quantiques.
Sinon, donne moi un exemple de méthode de tirage direct différent de la méthode rand, s'il te plait.

Est-ce vraiment nécessaire d'insister ?

Puisque tu le souhaites... Lecture du dernier bit du dernier registre visité : tirage binaire.

Et voilà.

Et soit dit en passant, le problème de gloups13 est réglé, quant à Bayes, comprenne qui pourra.

[minushabens]

en R, si le vecteur des probabilités est x, on obtient un tirage simplement par la commande

sum(cumsum(x)<runif(1))+1

évidemment si on doit répéter de nombreuses fois ce tirage on a intérêt à calculer la fonction de répartition (cumsum) une fois pour toutes.

[invite73192618]

Al là je n'y arrive pas. Car si P_A=10^-8 ,comment faire un programme qui tienne compte de ce cas. Je vois pas comment faire.

Ton probleme est regle plus haut, mais je pense qu'il y a une facon de faire possiblement plus proche de ta logique initiale.

1) on tire au hasard dans quel ordre on examine les possibilites
2) on teste si la probabilite est superieur a un critere initialise a 1
3a) si non on diminue le critere de la probabilite associee a la possibilite examinee et on examine la suivante
3b) si oui la possibilite courante est retenue

ordre=randperm(190); qui=1; critere=1;
While b(ordre(qui))<critere
        qui=qui+1; critere=critere-b(ordre(qui));
end 
disp(num2str(ordre(qui))

Il pourrait y avoir des cas dans lesquels cette methode est plus rapide, en particulier si b prend beaucoup de place en memoire pour de mauvaises raisons (exemple: si la plupart des probabilites ont deux ou trois chiffres de precision et une seule demande 8 chiffres, alors la methode vectorielle impose que toutes les probabilites soient decrites avec un trop grand nombre de chiffres), ou si b est tellement grand qu'il ne tient pas en memoire vive, auquel cas cette derniere methode evite (en moyenne) de manipuler la moitier des elements de b.

Une amelioration supplementaire serait de generer le vecteur 'ordre' a la volee, mais je ne sais pas s'il y a une methode elegante pour faire ca sous matlab. Evil?

[Evil.Saien]

Humhum...

Juste pour voir, teste ton algorithme avec 2 boules, la boule A avec pA = 0.1 et la boule B avec pB = 0.9 (donc b = [0.1 0.9])

ordre est soit [1 2] soit [2 1]

Si c'est [1 2] le while stoppe à la 2ème itération et te retourne 2
Si c'est [2 1] le while stoppe à la 2ème itération également et te retourne 1

Donc t'as perdu la pondération du tirage.

[invite73192618]

Tester un code avant de le poster? Quelle drôle d'idée!

Blague à part je l'ai testé et...

nb=2; b=rand(nb,1); b=b/sum(b);
ordre=randperm(nb); qui=1; critere=1;
while b(ordre(qui))>critere
        qui=qui+1; critere=critere-b(ordre(qui));
end 
disp(num2str(ordre(qui)))

...à part une majuscule malencontreuse au While, je ne retrouve pas l'erreur que tu signales.

[invite73192618]

Ah oui, pb avec le while, dsl je reviens.

[invite73192618]

...et avec le rand ...et avec le qui+1 ...et avec le break ...voilà une version apparemment ok, en mode verbeux pour bien voir ce qui se passe.

nb=190;
b=rand(nb,1); b=b/sum(b);
ordre=randperm(nb); qui=1; critere=rand(1); 
disp(['Étape #' num2str(qui) ' - ID de b #' num2str(ordre(qui)) ' - critere:' num2str(critere) ' - valeur de b: ' num2str(b(ordre(qui))) ])
while b(ordre(qui))<critere
        critere=critere-b(ordre(qui)); qui=qui+1; 
        disp(['Étape #' num2str(qui) ' - ID de b #' num2str(ordre(qui)) ' - critere:' num2str(critere) ' - valeur de b: ' num2str(b(ordre(qui))) ])
        if qui==nb; break; end
end 
disp(['Sélection du b #' num2str(ordre(qui))])

Étape #1 - ID de b #105 - critere:0.10436- valeur de b: 0.0011286
Étape #2 - ID de b #80 - critere:0.10323- valeur de b: 0.01031
Étape #3 - ID de b #145 - critere:0.09292- valeur de b: 0.0091743
Étape #4 - ID de b #98 - critere:0.083746- valeur de b: 0.005859
Étape #5 - ID de b #114 - critere:0.077887- valeur de b: 0.0083457
Étape #6 - ID de b #131 - critere:0.069541- valeur de b: 0.0065823
Étape #7 - ID de b #146 - critere:0.062959- valeur de b: 0.0025075
Étape #8 - ID de b #181 - critere:0.060451- valeur de b: 0.0068357
Étape #9 - ID de b #173 - critere:0.053616- valeur de b: 0.0055212
Étape #10 - ID de b #41 - critere:0.048095- valeur de b: 0.010026
Étape #11 - ID de b #188 - critere:0.038068- valeur de b: 0.007844
Étape #12 - ID de b #45 - critere:0.030224- valeur de b: 0.00015499
Étape #13 - ID de b #182 - critere:0.030069- valeur de b: 0.0028805
Étape #14 - ID de b #180 - critere:0.027189- valeur de b: 0.0073223
Étape #15 - ID de b #108 - critere:0.019866- valeur de b: 0.0056085
Étape #16 - ID de b #16 - critere:0.014258- valeur de b: 0.0033861
Étape #17 - ID de b #158 - critere:0.010872- valeur de b: 0.0044176
Étape #18 - ID de b #58 - critere:0.006454- valeur de b: 0.0058153
Étape #19 - ID de b #160 - critere:0.00063869- valeur de b: 0.00090526
Sélection du b #160
>>

Donc retour à ma question: une méthode élégante pour générer le vecteur 'ordre' a la volée plutôt qu'en initialisation de programme?

[Evil.Saien]

Ok, ça marche mais c'est uniquement dû au rand introduit

Le randperm ne sert juste... à rien.

Voici comment le vérifier :

nb=3;
b=[0.1 0.1 0.8]; b=b/sum(b);
tirages = 100000;
ordre = 1:nb;
resultat = zeros(1, tirages);
for ii = 1 : tirages
    qui=1; critere=rand(1); 
    while b(ordre(qui))<critere
        critere=critere-b(ordre(qui)); qui=qui+1; 
        if qui==nb; break; end
    end
    resultat(ii) = ordre(qui);
end

Il faut fait beaucoup de tirages pour être sûr que les pondérations sont gardées. Ensuite la commande

tabulate(resultat)

te donne les probabilités sorties:

>> tabulate(resultat)
  Value    Count   Percent
      1     9995      9.99%
      2     9917      9.92%
      3    80088     80.09%

Bref, c'est exactement

find(cumsum(b)<rand(1)), 1, 'last')

[invite73192618]

Bref, c'est exactement

find(cumsum(b)<rand(1)), 1, 'last')

Très joli, mais non ce n'est pas vraiment la même chose. Le résultat est bien le même (heureusement!), mais la complexité de l'algorithme est différente. Quand tu fais cumsum (et peut-être find aussi mais je ne suis pas certain), tu imposes de lire la totalité de b, et donc la vitesse est au minimum fonction linéaire de la taille de b. Au contraire l'algorithme précédent n'est pas juste beaucoup plus laid à écrire, il est aussi plus performant quand la taille de b devient inconfortablement grande.

Règle sur le pouce: tant que b rentre en mémoire vive, l'écriture typique matlab que tu indiques sera la meilleure. Mais dès que b est trop gros il faut devenir très prudent avec l'écriture typique matlab. Cela dit, quand on manipule des trucs plus gros que la mémoire vive, une autre règle sur le pouce c'est de passer à python.

Le randperm ne sert juste... à rien.

Idem: randperm ne sert à rien au sens que le résultat est bien le même avec ou sans randperm, mais il permet de régulariser le comportement de l'algorithme quand la fonction de répartition est différente d'une loi uniforme (si tu retournes au message initial gloups13 mentionne des 10-2 et des 10-8, ce qui suggère une loi de puissance). Cela dit, ce n'est pas la meilleure façon de faire. Le meilleur serait d'échantillonner la base de donnée contenant les b(i) avec une marche aléatoire vers une cible fixée par rand(1). Pour que cela marche il faut absolument que celle-ci soit structurée (i.e. en ordre de probabilité croissante ou décroissante), mais avec cette condition je suis à peu près certain qu'on trouvera un algorithme dont le temps d'exécution varie en fonction du log de la taille de b (ce serait trivial si b contenait la probabilité cumulative mais je pense qu'on peut s'en sortir en traitant les probabilités individuelles directement).

[Evil.Saien]

Très joli, mais non ce n'est pas vraiment la même chose. Le résultat est bien le même (heureusement!), mais la complexité de l'algorithme est différente. Quand tu fais cumsum (et peut-être find aussi mais je ne suis pas certain), tu imposes de lire la totalité de b, et donc la vitesse est au minimum fonction linéaire de la taille de b. Au contraire l'algorithme précédent n'est pas juste beaucoup plus laid à écrire, il est aussi plus performant quand la taille de b devient inconfortablement grande

Faire un scan sur un vecteur, c'est une complexité linéaire donc c'est hyper rapide.

Mais bon, on peut tester avec un grand nb, beaucoup de tirages et voir ce qui est le plus rapide...

nb=2000;
b = rand(nb,1)'; b=b/sum(b);
tirages = 1000000;
resultat = zeros(1, tirages);
tic;
for ii = 1 : tirages
    ordre = randperm(nb); qui=1; critere=rand(1); 
    while b(ordre(qui))<critere
        critere=critere-b(ordre(qui)); qui=qui+1; 
        if qui==nb; break; end
    end
    resultat(ii) = ordre(qui);
end
toc

tic; resultat2 = sum(bsxfun(@gt, rand(tirages, 1), [0 cumsum(b)]), 2); toc

J'obtiens le résultat suivant:

Elapsed time is 200.072958 seconds.
Elapsed time is 6.850775 seconds.

Note que si tu n'utilises pas le randperm et définis b = 1:nb, ta méthode descend à 35 secondes environ, démontrant qu'un randperm à chaque itération est super couteux.