Je suis confronté au problème de programmation suivant. J'ai besoin de générer des n (a, b) tuples pour lesquels la somme de tous les a est une donnée A et la somme de tous les b est une donnée B et pour chaque tuple, le rapport de a / b est dans la plage (c_min, c_max). A / B se situe également dans la même plage. J'essaie également de m'assurer qu'il n'y a pas de biais dans le résultat autre que ce qui est introduit par les contraintes et que les valeurs a / b sont plus ou moins uniformément réparties dans la plage donnée.

Quelques clarifications et méta-contraintes:

  • A, B, c_min et c_max sont donnés.
  • Le rapport A / B se situe dans la plage (c_min, c_max). Il doit en être ainsi si le problème doit avoir une solution compte tenu des autres contraintes.
  • a et b sont >0 et non entiers.

J'essaie d'implémenter cela en Python mais les idées dans n'importe quelle langue (anglais inclus) sont très appréciées.

12
ktdrv 27 oct. 2011 à 00:52

5 réponses

Meilleure réponse

Commencez par générer autant de tuples identiques, n, que vous en avez besoin:

(A/n, B/n)

Maintenant, choisissez deux tuples au hasard. Apportez une modification aléatoire à la valeur a de l'un et une modification compensatoire à la valeur a de l'autre, tout en respectant les contraintes données. Remettez les deux tuples.

Maintenant, choisissez une autre paire aléatoire. Cette fois, tordez avec les valeurs b.

Faire mousser, rincer de nouveau.

2
rossum 26 oct. 2011 à 21:12

Je pense que la chose la plus simple est de

  1. Utilisez votre méthode préférée pour lancer des n-1 valeurs telles que \sum_i=0,n-1 a_i < A et définissez a_n pour obtenir le bon total. Il y a plusieurs questions à ce sujet, même si je n'ai jamais vu de réponse dont je suis vraiment satisfait. Je vais peut-être écrire un article ou quelque chose.

  2. Obtenez les n-1 b en lançant les c_i uniformément sur la plage autorisée, et définissez la finale b pour obtenir le bon total et vérifier le c final (je pense ça doit être OK, mais je ne l'ai pas encore prouvé).

Notez que puisque nous avons 2 contraintes dures, nous devons nous attendre à lancer 2n-2 nombres aléatoires, et cette méthode fait exactement cela (en supposant que vous pouvez faire l'étape 1 avec n-1 lancers.

2
dmckee --- ex-moderator kitten 26 oct. 2011 à 21:12

Voici donc ce que je pense du point de vue mathématique. Nous avons des séquences a_i et b_i telles que la somme de a_i est A et la somme de b_i est B. De plus, A/B est en (x,y), tout comme a_i/b_i pour chaque i. De plus, vous souhaitez que a_i/b_i soit uniformément distribué dans (x,y).

Faites-le donc à partir de la fin. Choisissez c_i dans (x,y) de sorte qu'ils soient uniformément distribués. Ensuite, nous voulons avoir l'égalité suivante a_i/b_i = c_i, donc a_i = b_i*c_i.

Il nous suffit donc de trouver b_i. Mais nous avons le système d'équations linéaires suivant:

A = (sum)b_i*c_i
B = (sum)b_i

b_i sont des variables. Résolvez-le (quelques astuces d'algèbre linéaire fantaisie) et vous avez terminé!

Notez que pour un assez grand n ce système aura beaucoup de solutions. Ils dépendront de certains paramètres que vous pouvez choisir au hasard.


Assez de l'approche théorique, voyons une solution pratique.

// EDIT 1: Voici du code Python noyau dur: D

import random
min = 0.0
max = 10.0
A = 500.0
B = 100.0

def generate(n):
    C = [min + i*(max-min)/(n+1) for i in range(1, n+1)]
    Y = [0]
    for i in range(1,n-1):
        # This line should be changed in order to always get positive numbers
        # It should be relatively easy to figure out some good random generator
        Y.append(random.random())
    val = A - C[0]*B
    for i in range(1, n-1):
        val -= Y[i] * (C[i] - C[0])
    val /= (C[n-1] - C[0])
    Y.append(val)
    val = B
    for i in range(1, n):
        val -= Y[i]
    Y[0] = val
    result = []
    for i in range(0, n):
        result.append([ Y[i]*C[i], Y[i] ])
    return result

Le résultat est une liste de paires (X,Y) répondant à vos conditions à l'exception qu'elles peuvent être négatives (voir la ligne du générateur aléatoire dans le code), c'est-à-dire que la première et la dernière paire peuvent contenir des nombres négatifs.

// EDIT 2:

Pour vous assurer qu'ils sont positifs, vous pouvez essayer quelque chose comme

Y.append(random.random() * B / n)

Au lieu de

Y.append(random.random())

Je ne suis pas sûr cependant.

// EDIT 3:

Afin d'avoir de meilleurs résultats, essayez quelque chose comme ceci:

avrg = B / n
ran = avrg / 20
for i in range(1, n-1):
    Y.append(random.gauss(avrg, ran))

Au lieu de

for i in range(1, n-1):
    Y.append(random.random())

Cela rendra tous les b_i proches de B / n. Malheureusement, le dernier mandat sautera encore parfois haut. Je suis désolé, mais il n'y a aucun moyen d'éviter cela (mathématiques) car le dernier et le premier termes dépendent des autres. Pour les petits n (~ 100), cela semble bien. Malheureusement, certaines valeurs négatives peuvent apparaître.

Le choix d'un générateur correct n'est pas si simple si vous souhaitez en outre que b_i soit uniformément distribué.

0
freakish 27 oct. 2011 à 11:17

L'échantillonnage de Gibbs bloqué est assez simple et converge vers la bonne distribution (c'est dans le sens de ce que propose Alexandre).

  1. Pour tout i, initialisez a i = A / n et b i = B / n.
  2. Sélectionnez i ≠ j uniformément au hasard. Avec la probabilité 1/2, mettez à jour un i et un j avec des valeurs aléatoires uniformes satisfaisant les contraintes. Le reste du temps, faites de même pour b i et b j .
  3. Répétez l'étape 2 autant de fois que cela semble nécessaire pour votre application. Je n'ai aucune idée du taux de convergence.
1
anononononononon 27 oct. 2011 à 14:25

Beaucoup de bonnes idées ici. Merci! L'idée de Rossum semblait la plus simple à mettre en œuvre, alors je me suis lancée. Voici le code de la postérité:

c_min = 0.25
c_max = 0.75
a_sum = 100.0
b_sum = 200.0
n = 1000 

a = [a_sum / n] * n
b = [b_sum / n] * n

while not good_enough(a, b):
    i, j = random.sample(range(n), 2)
    li, ui = c_min * b[i] - a[i], c_max * b[i] - a[i]
    lj, uj = a[j] - c_min * b[j], a[j] - c_max * b[j]
    llim = max((li, uj))
    ulim = min((ui, lj))
    q = random.uniform(llim, ulim)
    a[i] += q
    a[j] -= q

    i, j = random.sample(range(n), 2)
    li, ui = a[i] / c_max - b[i], a[i] / c_min - b[i]
    lj, uj = b[j] - a[j] / c_max, b[j] - a[j] / c_min
    llim = max((li, uj))
    ulim = min((ui, lj))
    q = random.uniform(llim, ulim)
    b[i] += q
    b[j] -= q

La fonction good_enough(a, b) peut être beaucoup de choses. J'ai essayé:

  • L'écart type, qui est aléatoire, car vous ne savez pas ce qu'est une valeur suffisamment bonne.
  • Kurtosis, où une grande valeur négative serait bien. Cependant, il est relativement lent à calculer et n'est pas défini avec les valeurs de départ de (a_sum / n, b_sum / n) (bien que ce soit trivial à corriger).
  • Asymétrie, où une valeur proche de 0 est souhaitable. Mais il présente les mêmes inconvénients que le kurtosis.
  • Un certain nombre d'itérations proportionnelles à n. 2n parfois n'était pas suffisant, n ^ 2 est un peu exagéré et, bien, exponentiel.

Idéalement, une heuristique utilisant une combinaison d'asymétrie et de kurtosis serait la meilleure, mais je me suis contenté de m'assurer que chaque valeur a été modifiée par rapport à l'initiale (encore une fois, comme rossum l'a suggéré dans un commentaire). Bien qu'il n'y ait aucune garantie théorique que la boucle se terminera, cela a semblé fonctionner assez bien pour moi.

0
ktdrv 28 oct. 2011 à 01:03