Générer des nombres aléatoires avec une distribution probabiliste

Question

Ok, alors voici mon problème. Nous cherchons à acheter un ensemble de données d'une entreprise pour augmenter notre ensemble de données existant. Pour les besoins de cette question, disons que cet ensemble de données classe les places avec un nombre organique (ce qui signifie que le nombre attribué à un endroit n'a aucune incidence sur le nombre attribué à un autre). La plage technique est de 0 à l'infini, mais à partir des échantillons que j'ai vus, il est de 0 à 70. Sur la base de l'échantillon, ce n'est certainement pas une distribution uniforme (sur 10 000, il y en a peut-être 5 50 avec un score supérieur à 10 et 1000 avec un score supérieur à 1). Avant que nous décidions d'acheter cet ensemble, nous aimerions le simuler afin que nous puissions voir à quel point il peut être utile. Donc, pour le simuler, j'ai pensé à générer un nombre aléatoire pour chaque lieu (environ 150 000 nombres aléatoires). Mais, je veux aussi garder l'esprit des données, et garder la distribution relativement la même (ou du moins raisonnablement proche). Je me suis creusé la tête toute la journée en essayant de trouver un moyen de le faire, et je suis arrivé vide. Une idée que j'avais pensée était de mettre au carré le nombre aléatoire (entre 0 et sqrt (70)). Mais cela favoriserait à la fois moins de 1 et plus grands nombres. Je pense que la distribution réelle devrait être hyperbolique dans le premier quadrant ... Je suis juste en train de masquer comment transformer une distribution linéaire et régulière de nombres aléatoires en une distribution hyperbolique (Si hyperbolique est même ce que je vouloir en premier lieu).

Des pensées?

Donc, pour résumer, voici la distribution je voudrais (environ):

• 40 - 70: 0,02% - 0,05%
• 10-40: 0,5% - 1%
• 1 - 10: 10% - 20%
• 0 - 1.: Remainder (78,95% - 89,48%)

Answer 1

Regardez les distributions utilisées dans l'analyse de fiabilité - elles ont tendance à avoir ces longues queues. Une possibilité relativement simple est la distribution de Weibull avec P (X> x) = exp [- (x/b)^a]. Si vous placez vos valeurs comme P (X> 1) = 0.1 et P (X> 10) = 0.005, j'ai = 0.36 et b = 0.1. Cela impliquerait que P (X> 40) * 10000 = 1,6, ce qui est un peu trop bas, mais P (X> 70) * 10000 = 0,2 ce qui est raisonnable.

EDIT Oh, et pour générer une variable aléatoire à distribution de Weibull à partir d'un uniforme (0,1) Valeur U, calculer simplement b * [- log (1-u)]^(1/a). C'est la fonction inverse de 1-P (X> x) au cas où j'ai mal calculé quelque chose.

Answer 2

il y a quelques années écrites pour PHP4, il suffit de choisir votre distribution:

<?php 

define('RandomGaussian',   'gaussian') ;   // gaussianWeightedRandom() 
define('RandomBell',    'bell') ;    // bellWeightedRandom() 
define('RandomGaussianRising',  'gaussianRising') ; // gaussianWeightedRisingRandom() 
define('RandomGaussianFalling', 'gaussianFalling') ; // gaussianWeightedFallingRandom() 
define('RandomGamma',    'gamma') ;    // gammaWeightedRandom() 
define('RandomGammaQaD',   'gammaQaD') ;   // QaDgammaWeightedRandom() 
define('RandomLogarithmic10',  'log10') ;    // logarithmic10WeightedRandom() 
define('RandomLogarithmic',  'log') ;    // logarithmicWeightedRandom() 
define('RandomPoisson',   'poisson') ;   // poissonWeightedRandom() 
define('RandomDome',    'dome') ;    // domeWeightedRandom() 
define('RandomSaw',    'saw') ;    // sawWeightedRandom() 
define('RandomPyramid',   'pyramid') ;   // pyramidWeightedRandom() 
define('RandomLinear',    'linear') ;   // linearWeightedRandom() 
define('RandomUnweighted',   'non') ;    // nonWeightedRandom() 



function mkseed() 
{ 
    srand(hexdec(substr(md5(microtime()), -8)) & 0x7fffffff) ; 
} // function mkseed() 




/* 
function factorial($in) { 
    if ($in == 1) { 
     return $in ; 
    } 
    return ($in * factorial($in - 1.0)) ; 
} // function factorial() 


function factorial($in) { 
    $out = 1 ; 
    for ($i = 2; $i <= $in; $i++) { 
     $out *= $i ; 
    } 

    return $out ; 
} // function factorial() 
*/ 




function random_0_1() 
{ 
    // returns random number using mt_rand() with a flat distribution from 0 to 1 inclusive 
    // 
    return (float) mt_rand()/(float) mt_getrandmax() ; 
} // random_0_1() 


function random_PN() 
{ 
    // returns random number using mt_rand() with a flat distribution from -1 to 1 inclusive 
    // 
    return (2.0 * random_0_1()) - 1.0 ; 
} // function random_PN() 




function gauss() 
{ 
    static $useExists = false ; 
    static $useValue ; 

    if ($useExists) { 
     // Use value from a previous call to this function 
     // 
     $useExists = false ; 
     return $useValue ; 
    } else { 
     // Polar form of the Box-Muller transformation 
     // 
     $w = 2.0 ; 
     while (($w >= 1.0) || ($w == 0.0)) { 
      $x = random_PN() ; 
      $y = random_PN() ; 
      $w = ($x * $x) + ($y * $y) ; 
     } 
     $w = sqrt((-2.0 * log($w))/$w) ; 

     // Set value for next call to this function 
     // 
     $useValue = $y * $w ; 
     $useExists = true ; 

     return $x * $w ; 
    } 
} // function gauss() 


function gauss_ms($mean, 
        $stddev) 
{ 
    // Adjust our gaussian random to fit the mean and standard deviation 
    // The division by 4 is an arbitrary value to help fit the distribution 
    //  within our required range, and gives a best fit for $stddev = 1.0 
    // 
    return gauss() * ($stddev/4) + $mean; 
} // function gauss_ms() 


function gaussianWeightedRandom($LowValue, 
           $maxRand, 
           $mean=0.0, 
           $stddev=2.0) 
{ 
    // Adjust a gaussian random value to fit within our specified range 
    //  by 'trimming' the extreme values as the distribution curve 
    //  approaches +/- infinity 
    $rand_val = $LowValue + $maxRand ; 
    while (($rand_val < $LowValue) || ($rand_val >= ($LowValue + $maxRand))) { 
     $rand_val = floor(gauss_ms($mean,$stddev) * $maxRand) + $LowValue ; 
     $rand_val = ($rand_val + $maxRand)/2 ; 
    } 

    return $rand_val ; 
} // function gaussianWeightedRandom() 


function bellWeightedRandom($LowValue, 
          $maxRand) 
{ 
    return gaussianWeightedRandom($LowValue, $maxRand, 0.0, 1.0) ; 
} // function bellWeightedRandom() 


function gaussianWeightedRisingRandom($LowValue, 
             $maxRand) 
{ 
    // Adjust a gaussian random value to fit within our specified range 
    //  by 'trimming' the extreme values as the distribution curve 
    //  approaches +/- infinity 
    // The division by 4 is an arbitrary value to help fit the distribution 
    //  within our required range 
    $rand_val = $LowValue + $maxRand ; 
    while (($rand_val < $LowValue) || ($rand_val >= ($LowValue + $maxRand))) { 
     $rand_val = $maxRand - round((abs(gauss())/4) * $maxRand) + $LowValue ; 
    } 

    return $rand_val ; 
} // function gaussianWeightedRisingRandom() 


function gaussianWeightedFallingRandom($LowValue, 
             $maxRand) 
{ 
    // Adjust a gaussian random value to fit within our specified range 
    //  by 'trimming' the extreme values as the distribution curve 
    //  approaches +/- infinity 
    // The division by 4 is an arbitrary value to help fit the distribution 
    //  within our required range 
    $rand_val = $LowValue + $maxRand ; 
    while (($rand_val < $LowValue) || ($rand_val >= ($LowValue + $maxRand))) { 
     $rand_val = floor((abs(gauss())/4) * $maxRand) + $LowValue ; 
    } 

    return $rand_val ; 
} // function gaussianWeightedFallingRandom() 


function logarithmic($mean=1.0, $lambda=5.0) 
{ 
    return ($mean * -log(random_0_1()))/$lambda ; 
} // function logarithmic() 


function logarithmicWeightedRandom($LowValue, 
            $maxRand) 
{ 
    do { 
     $rand_val = logarithmic() ; 
    } while ($rand_val > 1) ; 

    return floor($rand_val * $maxRand) + $LowValue ; 
} // function logarithmicWeightedRandom() 


function logarithmic10($lambda=0.5) 
{ 
    return abs(-log10(random_0_1())/$lambda) ; 
} // function logarithmic10() 


function logarithmic10WeightedRandom($LowValue, 
             $maxRand) 
{ 
    do { 
     $rand_val = logarithmic10() ; 
    } while ($rand_val > 1) ; 

    return floor($rand_val * $maxRand) + $LowValue ; 
} // function logarithmic10WeightedRandom() 


function gamma($lambda=3.0) 
{ 
    $wLambda = $lambda + 1.0 ; 
    if ($lambda <= 8.0) { 
     // Use direct method, adding waiting times 
     $x = 1.0 ; 
     for ($j = 1; $j <= $wLambda; $j++) { 
      $x *= random_0_1() ; 
     } 
     $x = -log($x) ; 
    } else { 
     // Use rejection method 
     do { 
      do { 
       // Generate the tangent of a random angle, the equivalent of 
       //  $y = tan(pi * random_0_1()) 
       do { 
        $v1 = random_0_1() ; 
        $v2 = random_PN() ; 
       } while (($v1 * $v1 + $v2 * $v2) > 1.0) ; 
       $y = $v2/$v1 ; 
       $s = sqrt(2.0 * $lambda + 1.0) ; 
       $x = $s * $y + $lambda ; 
      // Reject in the region of zero probability 
      } while ($x <= 0.0) ; 
      // Ratio of probability function to comparison function 
      $e = (1.0 + $y * $y) * exp($lambda * log($x/$lambda) - $s * $y) ; 
     // Reject on the basis of a second uniform deviate 
     } while (random_0_1() > $e) ; 
    } 

    return $x ; 
} // function gamma() 


function gammaWeightedRandom($LowValue, 
           $maxRand) 
{ 
    do { 
     $rand_val = gamma()/12 ; 
    } while ($rand_val > 1) ; 

    return floor($rand_val * $maxRand) + $LowValue ; 
} // function gammaWeightedRandom() 


function QaDgammaWeightedRandom($LowValue, 
           $maxRand) 
{ 
    return round((asin(random_0_1()) + (asin(random_0_1()))) * $maxRand/pi()) + $LowValue ; 
} // function QaDgammaWeightedRandom() 


function gammaln($in) 
{ 
    $tmp = $in + 4.5 ; 
    $tmp -= ($in - 0.5) * log($tmp) ; 

    $ser = 1.000000000190015 
      + (76.18009172947146/$in) 
      - (86.50532032941677/($in + 1.0)) 
      + (24.01409824083091/($in + 2.0)) 
      - (1.231739572450155/($in + 3.0)) 
      + (0.1208650973866179e-2/($in + 4.0)) 
      - (0.5395239384953e-5/($in + 5.0)) ; 

    return (log(2.5066282746310005 * $ser) - $tmp) ; 
} // function gammaln() 


function poisson($lambda=1.0) 
{ 
    static $oldLambda ; 
    static $g, $sq, $alxm ; 

    if ($lambda <= 12.0) { 
     // Use direct method 
     if ($lambda <> $oldLambda) { 
      $oldLambda = $lambda ; 
      $g = exp(-$lambda) ; 
     } 
     $x = -1 ; 
     $t = 1.0 ; 
     do { 
      ++$x ; 
      $t *= random_0_1() ; 
     } while ($t > $g) ; 
    } else { 
     // Use rejection method 
     if ($lambda <> $oldLambda) { 
      $oldLambda = $lambda ; 
      $sq = sqrt(2.0 * $lambda) ; 
      $alxm = log($lambda) ; 
      $g = $lambda * $alxm - gammaln($lambda + 1.0) ; 
     } 
     do { 
      do { 
       // $y is a deviate from a Lorentzian comparison function 
       $y = tan(pi() * random_0_1()) ; 
       $x = $sq * $y + $lambda ; 
      // Reject if close to zero probability 
      } while ($x < 0.0) ; 
      $x = floor($x) ; 
      // Ratio of the desired distribution to the comparison function 
      // We accept or reject by comparing it to another uniform deviate 
      // The factor 0.9 is used so that $t never exceeds 1 
      $t = 0.9 * (1.0 + $y * $y) * exp($x * $alxm - gammaln($x + 1.0) - $g) ; 
     } while (random_0_1() > $t) ; 
    } 

    return $x ; 
} // function poisson() 


function poissonWeightedRandom($LowValue, 
           $maxRand) 
{ 
    do { 
     $rand_val = poisson()/$maxRand ; 
    } while ($rand_val > 1) ; 

    return floor($x * $maxRand) + $LowValue ; 
} // function poissonWeightedRandom() 


function binomial($lambda=6.0) 
{ 
} 


function domeWeightedRandom($LowValue, 
          $maxRand) 
{ 
    return floor(sin(random_0_1() * (pi()/2)) * $maxRand) + $LowValue ; 
} // function bellWeightedRandom() 


function sawWeightedRandom($LowValue, 
          $maxRand) 
{ 
    return floor((atan(random_0_1()) + atan(random_0_1())) * $maxRand/(pi()/2)) + $LowValue ; 
} // function sawWeightedRandom() 


function pyramidWeightedRandom($LowValue, 
           $maxRand) 
{ 
    return floor((random_0_1() + random_0_1())/2 * $maxRand) + $LowValue ; 
} // function pyramidWeightedRandom() 


function linearWeightedRandom($LowValue, 
           $maxRand) 
{ 
    return floor(random_0_1() * ($maxRand)) + $LowValue ; 
} // function linearWeightedRandom() 


function nonWeightedRandom($LowValue, 
          $maxRand) 
{ 
    return rand($LowValue,$maxRand+$LowValue-1) ; 
} // function nonWeightedRandom() 




function weightedRandom($Method, 
         $LowValue, 
         $maxRand) 
{ 
    switch($Method) { 
     case RandomGaussian   : 
      $rVal = gaussianWeightedRandom($LowValue, $maxRand) ; 
      break ; 
     case RandomBell    : 
      $rVal = bellWeightedRandom($LowValue, $maxRand) ; 
      break ; 
     case RandomGaussianRising : 
      $rVal = gaussianWeightedRisingRandom($LowValue, $maxRand) ; 
      break ; 
     case RandomGaussianFalling : 
      $rVal = gaussianWeightedFallingRandom($LowValue, $maxRand) ; 
      break ; 
     case RandomGamma   : 
      $rVal = gammaWeightedRandom($LowValue, $maxRand) ; 
      break ; 
     case RandomGammaQaD   : 
      $rVal = QaDgammaWeightedRandom($LowValue, $maxRand) ; 
      break ; 
     case RandomLogarithmic10 : 
      $rVal = logarithmic10WeightedRandom($LowValue, $maxRand) ; 
      break ; 
     case RandomLogarithmic  : 
      $rVal = logarithmicWeightedRandom($LowValue, $maxRand) ; 
      break ; 
     case RandomPoisson   : 
      $rVal = poissonWeightedRandom($LowValue, $maxRand) ; 
      break ; 
     case RandomDome    : 
      $rVal = domeWeightedRandom($LowValue, $maxRand) ; 
      break ; 
     case RandomSaw    : 
      $rVal = sawWeightedRandom($LowValue, $maxRand) ; 
      break ; 
     case RandomPyramid   : 
      $rVal = pyramidWeightedRandom($LowValue, $maxRand) ; 
      break ; 
     case RandomLinear   : 
      $rVal = linearWeightedRandom($LowValue, $maxRand) ; 
      break ; 
     default      : 
      $rVal = nonWeightedRandom($LowValue, $maxRand) ; 
      break ; 
    } 

    return $rVal; 
} 

?> 

Answer 3

Cette manière naïve de le faire va probablement fausser la distribution d'une manière que je ne peux pas voir maintenant. L'idée est simplement d'itérer sur votre premier ensemble de données, triées et en paires. Ensuite, randomize 15 nouveaux numéros entre chaque paire pour obtenir le nouveau tableau.

Ruby exemple, puisque je ne parle pas beaucoup de PHP. Espérons qu'une telle idée simple devrait être facile à traduire en PHP.

numbers=[0.1,0.1,0.12,0.13,0.15,0.17,0.3,0.4,0.42,0.6,1,3,5,7,13,19,27,42,69] 
more_numbers=[] 
numbers.each_cons(2) { |a,b| 15.times { more_numbers << a+rand()*(b-a) } } 
more_numbers.sort! 

Answer 4

Le plus facile (mais pas très efficace) façon de générer des nombres aléatoires qui suivent une distribution donnée est une technique appelée Von Neumann Rejection.

La simple explication de la technique est la suivante. Créez une boîte qui entoure complètement votre distribution. (Appelez votre distribution f) Ensuite, choisissez un point aléatoire (x,y) dans la boîte. Si y < f(x), utilisez x comme un nombre aléatoire. Si y > f(x), alors ignorez x et y et choisissez un autre point. Continuez jusqu'à ce que vous ayez une quantité suffisante de valeurs à utiliser. Les valeurs de x que vous ne rejetez pas seront distribuées selon f.