C – générer des nombres aléatoires dans un intervalle par rapport à une moyenne

Je dois générer un ensemble de nombres aléatoires dans un intervalle qui a également une valeur moyenne. Par exemple, min = 1000, max = 10000 et une moyenne de 7000. Je sais comment créer des nombres dans une plage, mais je me bats avec le facteur de valeur moyenne. Y a-t-il une fonction que je peux utiliser?

Ce que vous recherchez se fait le plus facilement avec la méthode dite de refus d’acceptation.

Divisez votre intervalle en intervalles plus petits. Spécifier une fonction de densité de probabilité (PDF) peut être très simple aussi, comme une fonction de pas. Pour la dissortingbution gaussienne, vous auriez des marches à gauche et à droite plus basses que votre marche du milieu, c.-à-d.

Approche générale pour générer des valeurs aléatoires sur un fichier PDF

Générez un nombre aléatoire dans tout l’intervalle. Si le nombre généré est supérieur à la valeur de votre PDF à ce stade, rejetez le nombre généré.

Répétez les étapes jusqu’à obtenir le nombre de points souhaité.


EDIT 1

Preuve de concept sur un PDF gaussien.

Ok, donc l’idée de base est montrée dans le graphique (a).

  1. Définissez / choisissez votre fonction de densité de probabilité (PDF). PDF est une fonction statistique d’une variable aléatoire et décrit la probabilité de trouver la valeur x dans une mesure / expérience. Une fonction peut être un fichier PDF d’une variable aléatoire x si elle satisfait: 1) f(x) >= 0 et 2) est normalisée (c’est-à-dire qu’elle sum ou intègre jusqu’à la valeur 1).
  2. Obtenir le maximum ( max ) et les “zéro points” ( z1 < z2 ) du fichier PDF. Certains PDF peuvent avoir leurs points zéro à l'infini. Dans ce cas, déterminez les points de coupure (z1, z2) pour lesquels PDF(z1>x>z2) < eta où vous choisissez vous-même eta . Fondamentalement, cela signifie que vous définissez une valeur petite eta puis que vos points zéro sont les valeurs pour lesquelles la valeur de PDF(x) est inférieure à celle de eta.
  3. Définissez l'intervalle Ch(z1, z2, max) de votre générateur aléatoire. C'est l'intervalle dans lequel vous générez vos variables aléatoires.
  4. Générez une variable aléatoire x telle que z1 .
  5. Générez une seconde variable aléatoire y sans rapport dans la plage (0, max) . Si la valeur de y est inférieure à PDF(x) rejetez les deux valeurs générées de manière aléatoire (x,y) et revenez à l'étape 4. Si la valeur générée y est supérieure à PDF(x) acceptez la valeur x tant que point généré de manière aléatoire. sur une dissortingbution et le return .

Voici le code qui reproduit un comportement similaire pour un PDF gaussien.

 #include "Random.h" #include  using namespace std; double gaus(double a, double b, double c, double x) { return a*exp( -((xb)*(xb)/(2*c*c) )); } double* random_on_a_gaus_dissortingbution(double inter_a, double inter_b) { double res [2]; double a = 1.0; //currently parameters for the Gaussian double b = 2.0; //are defined here to avoid having double c = 3.0; //a long function declaration line. double x = kiss::Ran(inter_a, inter_b); double y = kiss::Ran(0.0, 1.0); while (y>gaus(a,b,c,x)) //keep creating values until step 5. is satisfied. { x = kiss::Ran(inter_a, inter_b); //this is interval (z1, z2) y = kiss::Ran(0.0, 1.0); //this is the interval (0, max) } res[0] = x; res[1] = y; return res; //I return (x,y) for plot reasons, only x is the randomly } //generated value you're looking for. void main() { double* x; ofstream f; f.open("test.txt"); for(int i=0; i<100000; i++) { //see bellow how I got -5 and 10 to be my interval (z1, z2) x = random_on_a_gaus_distribution(-5.0, 10.0); f << x[0]<<","< 

Étape 1

Nous définissons donc d’abord un aspect général d’un fichier PDF gaussien dans une fonction appelée gaus . Simple.

Ensuite, nous définissons une fonction random_on_a_gaus_dissortingbution qui utilise une fonction gaussienne bien définie. Dans une expérience \ mesure, nous obtiendrions les coefficients a, b, c en ajustant notre fonction. J'ai choisi des valeurs aléatoires (1, 2, 3) dans cet exemple. Vous pouvez choisir celles qui satisfont votre affectation matérielle (c'est-à-dire: les coefficients qui font une gaussienne dont la moyenne est de 7 000).

Étape 2 et 3

J'ai utilisé wolfram mathematica pour tracer gaus. avec les parameters 1,2,3 aussi voir quelles seraient les valeurs les plus appropriées pour max et (z1, z2) . Vous pouvez voir le graphique vous-même . Le maximum de la fonction est 1,0 et, selon une ancienne méthode scientifique appelée eyeballin, j’estimais que les seuils étaient compris entre -5,0 et 10,0.

Pour rendre random_on_a_gaus_dissortingbution plus générale, vous pouvez suivre l’étape 2) de manière plus rigoureuse et définir eta , puis calculer votre fonction en points successifs jusqu’à ce que PDF devienne plus petit que eta. Dangers avec ceci est que vos points de coupure peuvent être très éloignés et cela pourrait prendre longtemps pour des fonctions très monotones. De plus, vous devez trouver le maximum vous-même. Ceci est généralement délicat. Cependant, un problème plus simple est la minimisation du négatif d'une fonction. Cela peut aussi être délicat pour un cas général, mais pas "annulable". Le moyen le plus simple est de sortingcher un peu comme je l’ai fait et de coder ceci pour quelques fonctions seulement.

Étape 4 et 5

Et puis tu t'en vas. Continuez simplement à créer de nouveaux points jusqu'à ce que vous obteniez un résultat satisfaisant. VEUILLEZ NOTER que le nombre retourné x est un nombre aléatoire. Vous ne seriez pas en mesure de trouver un lien logique entre deux valeurs x créées successivement, ou le premier créé x et le millionième.

Cependant, le nombre de valeurs x acceptées dans l'intervalle autour de x_max de notre dissortingbution est supérieur au nombre de valeurs x créées dans des intervalles pour lesquels PDF(x) < PDF(x_max) .

Cela signifie simplement que vos nombres aléatoires seront pondérés dans l'intervalle choisi de manière à ce que la plus grande valeur PDF d'une variable aléatoire x corresponde à plus de points aléatoires acceptés dans un petit intervalle autour de cette valeur que autour de toute autre valeur de xi pour laquelle PDF(xi) .

Je suis retourné à la fois x et y pour pouvoir tracer le graphique ci-dessous, mais ce que vous cherchez à retourner est en réalité juste le x . J'ai fait les plots avec matplotlib.

Nuage de points de (x, y) valeurs, (random, probabilité_it_got_accepted_with)

Il est probablement préférable de n’afficher qu’un histogramme de variable créée aléatoirement sur une dissortingbution. Cela montre que les valeurs x proches de la valeur moyenne de votre fonction PDF sont les plus susceptibles d'être acceptées. Par conséquent, des variables créées de manière plus aléatoire avec ces valeurs approximatives seront créées.

Histogramme de la variable <code/> x </ code> juste créée de manière aléatoire dans la fonction <code> random_on_a_gaus_dissortingbution </ code>.

De plus, je suppose que vous seriez intéressé par la mise en oeuvre du générateur de nombre aléatoire kiss. IL EST TRÈS IMPORTANT D'AVOIR UN TRÈS BON GÉNÉRATEUR . J'ose dire que, dans une certaine mesure, kiss ne le coupe probablement pas (le mersene twister est souvent utilisé).

Random.h

 #pragma once #include  const unsigned RNG_MAX=4294967295; namespace kiss{ // unsigned int kiss_z, kiss_w, kiss_jsr, kiss_jcong; unsigned int RanUns(); void RunGen(); double Ran0(int upper_border); double Ran(double bottom_border, double upper_border); } namespace Crand{ double Ran0(int upper_border); double Ran(double bottom_border, double upper_border); } 

Kiss.cpp

 #include "Random.h" unsigned int kiss_z = 123456789; //od 1 do milijardu unsigned int kiss_w = 378295763; //od 1 do milijardu unsigned int kiss_jsr = 294827495; //od 1 do RNG_MAX unsigned int kiss_jcong = 495749385; //od 0 do RNG_MAX //KISS99* //Autor: George Marsaglia unsigned int kiss::RanUns() { kiss_z=36969*(kiss_z&65535)+(kiss_z>>16); kiss_w=18000*(kiss_w&65535)+(kiss_w>>16); kiss_jsr^=(kiss_jsr<<13); kiss_jsr^=(kiss_jsr>>17); kiss_jsr^=(kiss_jsr<<5); kiss_jcong=69069*kiss_jcong+1234567; return (((kiss_z<<16)+kiss_w)^kiss_jcong)+kiss_jsr; } void kiss::RunGen() { for (int i=0; i<2000; i++) kiss::RanUns(); } double kiss::Ran0(int upper_border) { unsigned velicinaIntervala = RNG_MAX / upper_border; unsigned granicaIzbora= velicinaIntervala*upper_border; unsigned slucajniBroj = kiss::RanUns(); while(slucajniBroj>=granicaIzbora) slucajniBroj = kiss::RanUns(); return slucajniBroj/velicinaIntervala; } double kiss::Ran (double bottom_border, double upper_border) { return bottom_border+(upper_border-bottom_border)*kiss::Ran0(100000)/(100001.0); } 

De plus, il existe les générateurs aléatoires standard C: CRands.cpp

 #include "Random.h" //standardni pseudo random generatori iz Ca double Crand::Ran0(int upper_border) { return rand()%upper_border; } double Crand::Ran (double bottom_border, double upper_border) { return (upper_border-bottom_border)*rand()/((double)RAND_MAX+1); } 

Il convient également de commenter le graphique (b) ci-dessus. Lorsque vous avez un PDF très mal comporté, PDF(x) variera considérablement entre les grands et les très petits nombres.

Le problème, c'est que la zone d'intervalle Ch(x) correspond bien aux valeurs extrêmes du PDF, mais puisque nous créons également une variable aléatoire y pour les petites valeurs de PDF(x) ; les chances d'accepter cette valeur sont infimes! Il est plus probable que la valeur y générée sera toujours plus grande que PDF(x) à ce stade. Cela signifie que vous passerez beaucoup de cycles à créer des nombres qui ne seront pas choisis et que tous vos nombres aléatoires choisis seront très localement liés au max de votre PDF.

C'est pourquoi il est souvent utile de ne pas avoir les mêmes intervalles Ch(x) partout, mais de définir un ensemble d'intervalles paramétrés. Cependant, cela ajoute un peu de complexité au code.

Où fixez-vous vos limites? Comment traiter les cas limites? Quand et comment déterminer que vous devez en effet soudainement utiliser cette approche? Le calcul de max peut ne pas être aussi simple à présent, cela dépend de la méthode que vous avez initialement envisagée.

De plus, vous devez maintenant corriger le fait que beaucoup plus de numéros sont acceptés plus facilement dans les zones où la hauteur de votre case Ch(x) est plus basse, ce qui biaise le PDF d'origine.

Cela peut être corrigé en pesant les nombres créés dans la limite inférieure par le rapport des hauteurs des limites supérieure et inférieure. En gros, répétez l'étape y une fois de plus. Créez un nombre aléatoire z compris entre 0 et 1 et comparez-le au ratio lower_height / higher_height, garanti <1. Si z est plus petit que le ratio: acceptez x et si c'est plus grand, rejetez.

Les généralisations du code présenté sont également possibles en écrivant une fonction qui prend à la place un pointeur d'object. En définissant votre propre classe, c’est-à-dire une function décrivant généralement les fonctions, ayant une méthode eval en un point, pouvant stocker vos parameters, calculer et stocker ses propres valeurs max / min et ses points zéro / limite, vous n’aurez pas à passer , ou les définir dans une fonction comme je l’ai fait.

Bonne chance, amuse toi bien!

tl; dr : élever une dissortingbution uniforme de 0 à 1 à la puissance (1 - m) / mm est la moyenne souhaitée (entre 0 et 1). Déplacer / mettre à l’échelle comme vous le souhaitez.


J’étais curieux de savoir comment mettre cela en œuvre. J’ai pensé qu’un trapèze serait la méthode la plus simple, mais vous êtes limité en ce sens que le moyen le plus extrême que vous puissiez obtenir est un sortingangle, qui n’est pas si extrême. Les maths ont commencé à devenir difficiles, alors je suis revenu à une méthode purement empirique qui semble bien fonctionner.

Quoi qu’il en soit, pour une dissortingbution, pourquoi ne pas commencer avec la dissortingbution uniforme [0, 1) et élever les valeurs à un pouvoir arbitraire. Place-les et la dissortingbution se déplace vers la droite. Racine carrée eux et ils décalent vers la gauche. Vous pouvez aller à l’extrême que vous voulez et pousser la dissortingbution aussi fort que vous le souhaitez.

 def randompow(p): return random.random() ** p 

(Tout est écrit en Python, mais devrait être assez facile à traduire. Si quelque chose n’est pas clair, il suffit de demander. random.random() renvoie les flottants de 0 à 1)

Alors, comment pouvons-nous ajuster ce pouvoir? Eh bien, comment la moyenne semble-t-elle évoluer avec des pouvoirs différents?

On dirait une sorte de courbe sigmoïde. Il y a beaucoup de fonctions sigmoïdes , mais la tangente hyperbolique semble plutôt bien fonctionner.

Pas à 100% là-bas, essayons de le redimensionner dans la direction X …

 # x are the values from -3 to 3 (log transformed from the powers used) # y are the empirically-determined means given all those powers def fitter(tanscale): xsc = tanscale * x sigtan = np.tanh(xsc) sigtan = (1 - sigtan) / 2 resid = sigtan - y return sum(resid**2) fit = scipy.optimize.minimize(fitter, 1) 

L’installateur indique que le meilleur facteur d’échelle est 1.1514088816214016. Les résidus sont en fait assez bas, alors ça sonne bien.

Implémenter l’inverse de toutes les maths dont je n’ai pas parlé ressemble à ceci:

 def distpow(mean): p = 1 - (mean * 2) p = np.arctanh(p) / 1.1514088816214016 return 10**p 

Cela nous donne le pouvoir d’utiliser dans la première fonction pour obtenir n’importe quel moyen de la dissortingbution. Une fonction usine peut renvoyer une méthode permettant de générer une série de nombres de la dissortingbution avec la moyenne souhaitée.

 def randommean(mean): p = distpow(mean) def f(): return random.random() ** p return f 

Comment ça va? Assez bien à 3-4 décimales:

 for x in [0.01, 0.1, 0.2, 0.4, 0.5, 0.6, 0.8, 0.9, 0.99]: f = randommean(x) # sample the dissortingbution 10 million times mean = np.mean([f() for _ in range(10000000)]) print('Target mean: {:0.6f}, actual: {:0.6f}'.format(x, mean)) Target mean: 0.010000, actual: 0.010030 Target mean: 0.100000, actual: 0.100122 Target mean: 0.200000, actual: 0.199990 Target mean: 0.400000, actual: 0.400051 Target mean: 0.500000, actual: 0.499905 Target mean: 0.600000, actual: 0.599997 Target mean: 0.800000, actual: 0.799999 Target mean: 0.900000, actual: 0.899972 Target mean: 0.990000, actual: 0.989996 

Une fonction plus succincte qui vous donne juste une valeur donnée à une moyenne (pas une fonction d’usine):

 def randommean(m): p = np.arctanh(1 - (2 * m)) / 1.1514088816214016 return random.random() ** (10 ** p) 

Edit: l’ ajustement avec le logarithme naturel de la moyenne au lieu de log10 a donné une valeur résiduelle suspecte proche de 0,5 Faire des calculs pour simplifier l’arctanh donne:

 def randommean(m): '''Return a value from the dissortingbution 0 to 1 with average *m*''' return random.random() ** ((1 - m) / m) 

À partir de là, il devrait être assez facile de changer, de redimensionner et d’arrondir la dissortingbution. La troncature en entier peut finir par décaler la moyenne de 1 (ou une demi-unité?), Donc c’est un problème non résolu (si c’est important).