12
votes

Comment générer des variables binaires corrélées

Je dois générer une série de variables binaires aléatoires n avec une fonction de corrélation donnée. Laissez x = { x i } est une série de variables binaires (prenant la valeur 0 ou 1, < em> i courir de 1 à n ). La probabilité marginale est donnée PR ( x i = 1) = p et les variables doivent être corrélées dans la manière suivante:

cor [ x i x j ] = const × | i - j | (pour i! = j)

où α est un nombre positif.

S'il est plus facile, considérez la fonction de corrélation:

cor [ x i x j ] = (| i - j | +1)

La partie essentielle est que je veux enquêter sur le comportement lorsque la fonction de corrélation va comme une loi de puissance. (pas α | i - j | )

est-il possible de générer une série comme celle-ci, de préférence en python?


13 commentaires

@Paul: Non, sans aucun doute pas: "L'objectif principal de Mathoverflow est que les utilisateurs de poser et de répondre aux questions de mathématiques de niveau de recherche, les types de questions que vous rencontrez lorsque vous écrivez ou lisez des articles ou des livres de niveau supérieurs." N'EST PAS TOYZ Ceci.


@Johannes: OK - Merci pour la clarification - cela semble être plus d'une question de mathématiques qu'une question de programmation en ce que la programmation sera triviale une fois que les mathématiques sont cloués, mais cela semble être plus approprié pour demander ce genre de question.


Est-il correct d'interpréter la corrélation comme à peu près "plus la différence des indices de deux variables est la plus petite est la probabilité de 1"? Intuitivement, cela préfère des grappes.


@Svante: dit d'une autre manière: plus la différence des indices des variables, plus ils se comportent comme des variables indépendantes. Votre intuition a raison, cela préfère les clusters.


Je ne pense pas que cela soit possible comme écrit. Lorsque N est très important, le produit attendu de XI largement séparé doit être très proche de zéro et nous pouvons trouver un ensemble de 1 / p + 1 qui doit essentiellement être mutuellement exclusif malgré la probabilité p. Peut-être que vous vouliez dire e [xi xj] - p ^ 2 = ...?


De plus, quelle gamme de n pensez-vous et acceptez-vous une fonction de décomposition différente (par exemple exponentielle, car ce serait beaucoup plus facile)?


@algorithmist: votre droit correct [XI XI] = E [XI XJ] - P ^ 2, j'espère avoir une N dans l'ordre de 1000-10000.


Cette configuration semble avoir beaucoup en commun avec le modèle de circuit binaire ISing, bien qu'il faudra des efforts pour déterminer exactement ce que votre structure de corrélation supposée implique. Parce que chaque article est connecté, je doute qu'il y ait une stratégie de simulation unique. Vous voulez probablement mettre en œuvre un échantillonneur GIBBS. Travaillez la distribution conditionnelle de tout élément particulier, tenant les autres fixes, puis simulez de manière itérative à partir de celle-ci, à partir d'une position aléatoire.


Jonalm: Êtes-vous sûr que c'est | i-j | -alpha plutôt que alpha ^ | i-j | ? (Juste vérifiant les fautes de frappe évidentes)


@jonalm: Pourriez-vous modifier votre question avec la formule de corrélation correcte? Il est plus facile de rechercher la fonction.


Jonalm, avez-vous besoin de mettre en œuvre cette structure de corrélation spécifique ou simplement une structure de corrélation dont la progéniture diminue avec la distance (et est modifiable par un paramètre)?


@Jason S, @aniko et @Matthhiu M. La caractéristique essentielle de la fonction de corrélation est la loi sur la puissance. J'ai édité la question. J'espère que c'est clair. @Tristan, c'est vrai, c'est essentiellement le modèle ising au point critique. Merci pour l'impuche.


@jonalm: Ok, je l'ai eu. Merci de clarifier.


6 Réponses :


0
votes

exprimer la distribution x i comme combinaison linéaire de quelques distributions de base indépendantes f j : > X I = A I1 F 1 + A I2 F 2 +. .. . Laissez-nous contraindre f j pour être des variables indépendantes uniformément distribuées en 0..1 ou en {0,1} (discret). Exprimons maintenant tout ce que nous savons dans la forme matricielle: xxx

Vous devez désormais résoudre les deux équations: xxx

solution de la première L'équation correspond à la recherche de la racine carrée de la matrice 3R ou 2R. Voir par exemple http://fr.wikipedia.org/wiki/cholesky_factorisation et généralement http://fr.wikipedia.org/wiki/square_root_of_a_matrix . Quelque chose devrait également être fait sur le second :)

Je demande aux mathématiciens de me corriger, car je pourrais très bien avoir mélangé à a avec un ou fait quelque chose encore plus faux.

pour générer une valeur de x i comme mélange linéaire des distributions de base, utilisez un processus en deux étapes: 1) Utilisez un Variable aléatoire uniforme Pour choisir l'une des distributions de base, pondérées avec une probabilité correspondante, 2) génère un résultat à l'aide de la distribution de base choisie.


5 commentaires

Malheureusement, la transition continue -> discrète est souvent la partie la plus difficile. Par exemple, le problème de la recherche d'une matrice Hadamard obtient un lot plus facile si des entrées complexes sont autorisées. Je ne vois aucun moyen de discrétiser votre solution dans le cadre donné.


Pourquoi devrait-il être difficile? La solution dépend simplement de la distribution résultante étant un mélange linéaire de distributions de base, mais je ne vois pas comment cela dépend de leur continuité. Est-ce que des distributions discrètes ne peuvent pas être facilement mélangées de manière linéaire?


Dans ce cas, c'est le fait que vos distributions continues ne constituent pas des combinaisons convexes des essais Bernoulli.


Je suis désolé, je ne comprends pas bien. Je sais quelle est une combinaison convexe et que sont les essais Bernoulli, mais toujours: j'ai édité mon poste; Le processus décrit n'est-il pas dans le dernier paragraphe donnant un résultat correct? Si oui, pourriez-vous me signaler à certaines sources élargir votre point? (Quoi qu'il en soit, je devrais probablement implémenter les choses que j'ai décrites et voir si cela fonctionne)


Le problème est que A a peut-être des entrées qui ne sont pas comprises entre 0 et 1.



0
votes

La solution de force brute consiste à exprimer les contraintes du problème en tant que programme linéaire avec 2 ^ n variables pr (w) w gammes sur toutes les chaînes binaires de longueur n . Premièrement, la contrainte que pr est une distribution de probabilité: xxx p> second, la contrainte que l'attente de chaque variable soit p : xxx

troisième, les contraintes de covariance: xxx

C'est très lent, mais une recherche de littérature superficielle n'a rien de mieux. Si vous décidez de la mettre en œuvre, voici quelques solvants LP avec des liaisons Python: http: // wiki. python.org/moin/numericandscientific/libraries


2 commentaires

Je ne connais pas la programmation linéaire, mais je ne vois pas comment cela fonctionnera. Toute configuration de la série binaire aura une probabilité non nulle. Est-il possible de calculer la probabilité de toute configuration?


Oui, si le problème est résoluble, la routine LP vous donnera la probabilité de chacune des configurations 2 ^ n.



1
votes

Une recherche rapide sur Resek révèle que r a des packages


2 commentaires

C'était aussi ma première idée, mais je doute qu'ils puissent gérer N = 1000-10000 002 par déclaration de Jonalm dans les commentaires.


Je ne sais pas r, mais cela semble être une raison d'apprendre :) Merci.



2
votes

Vous décrivez un processus aléatoire , et cela ressemble à un dur pour moi ... Si vous éliminez l'exigence binaire (0,1) et spécifiait la valeur et la variance attendues Il serait possible de décrire cela comme un générateur de bruit blanc d'alimentation à travers un filtre passe-bas 1 pôle, que je pense vous donnerait le α | ij | caractéristique.

Ceci pourrait réellement rencontrer la barre de Mathoverflow.net, en fonction de la façon dont il est formulé. Laissez-moi essayer de demander ....


mise à jour: j'ai fait Demandez sur Mathoverflow.net pour le α < sup> | ij | cas. Mais peut-être y avoir des idées qui peuvent être adaptées à votre cas.


4 commentaires

C'est en fait | I-J | ^ - \ alpha; la solution pour \ alpha ^ | i-j | est dans la littérature.


hmmm ... | i-j | -alpha n'a aucune solution pour i = j. Sommes-nous sûrs que l'OP n'a pas mal d'état?


Il peut être cor [xi xj] = const x | i-j | -alpha pour i! = J, ou cor [xj xi] = (| i-j | +1) ^ alfa (ce qui est le plus facile). Quoi qu'il en soit, je ne prétends pas qu'ils sont égaux, mais qui ne s'intéressait que par le comportement de la queue (| I-J | >> 1) de sorte que cela ne devrait pas avoir d'importance.


Merci beaucoup Jason. Bien que c'était une corrélation différente, les solutions étaient vraiment intéressantes.



0
votes

Voici une approche intuitive / expérimentale qui semble fonctionner.

si B est un R.V binaire, m est la moyenne du r.v. binaire, c est la corrélation que vous voulez, rand () génère un u (0,1) r.v., et d est le R.V binaire corrélé Vous voulez:

d = si (rand ()

c'est-à-dire si un r.v. uniforme est inférieur à la corrélation souhaitée, D = b. Sinon d = un autre nombre binaire aléatoire.

J'ai couru cela 1000 fois pour une colonne de 2000 binaires R.V.S. avec m = .5 et c = .4 et c = .5 La signification de la corrélation était exactement comme spécifiée, la distribution semblait normale. Pour une corrélation de 0,4, la déviation STD de la corrélation était de 0,02.

Désolé - je ne peux pas prouver que cela fonctionne tout le temps, mais vous devez admettre que c'est facile.


1 commentaires

Relire la question: ce n'est pas la bonne structure de corrélation.



5
votes

Merci pour toutes vos intrants. J'ai trouvé une réponse à ma question dans le mignon petit article de Chul Gyu Park et al., Donc, au cas où tout le monde se trouve dans le même problème, recherchez:

"Une méthode simple pour générer des variables binaires corrélées" (jstor.org.stable / 2684925)

pour un simple algorithme. L'algorithme fonctionne si tous les éléments de la matrice de corrélation sont positifs et pour une distribution marginale générale PR (X_I) = P_I.

j


2 commentaires

Intéressant ... Voir aussi projecteuclid.org/...


Avez-vous mis en œuvre cela à Python? Avez-vous toujours le script? Avez-vous une chance que vous puissiez l'inclure s'il vous plaît?