Genera numeri casuali con distribuzione probabilistica

Question

Ok, ecco il mio problema. Stiamo cercando di acquistare un set di dati da un'azienda per aumentare il nostro set di dati esistente. Ai fini di questa domanda, diciamo che questo set di dati classifica i luoghi con un numero organico (il che significa che il numero assegnato a un posto non ha alcun rapporto con il numero assegnato a un altro). Il range tecnico va da 0 a infinito, ma dai sample set che ho visto è compreso tra 0 e 70. Sulla base del campione, non è sicuramente una distribuzione uniforme (su 10.000 ci sono forse 5 posti con un punteggio superiore a 40, 50 con un punteggio superiore a 10 e 1000 con un punteggio superiore a 1). Prima di decidere di acquistare questo set, vorremmo simularlo in modo che possiamo vedere quanto possa essere utile.

Quindi, per simularlo, ho pensato di generare un numero casuale per ogni posto (circa 150.000 numeri casuali). Ma, voglio anche mantenere lo spirito dei dati e mantenere la distribuzione relativamente uguale (o almeno ragionevolmente vicina). Mi sto tormentando il cervello tutto il giorno cercando di pensare a un modo per farlo, e sono venuto fuori vuoto.

Un pensiero che avevo era di quadrare il numero casuale (tra 0 e sqrt (70)). Ma ciò favorirebbe meno di 1 e più numeri.

Sto pensando che la distribuzione reale dovrebbe essere iperbolica nel primo quadrante ... Sto solo ignorando come trasformare una distribuzione lineare e uniforme di numeri casuali in una distribuzione iperbolica (Se iperbolico è anche ciò che io voglio in primo luogo).

Qualche idea?

Quindi, per riassumere, ecco la distribuzione vorrei (circa):

• 40 - 70: 0,02% - 0,05%
• 10 - 40: 0,5% - 1%
• 1 - 10: 10% - 20%
• 0-1: Resto (78.95% - 89.48%)

Answer 1

Consultare le distribuzioni utilizzate nell'analisi dell'affidabilità - tendono ad avere queste lunghe code. Una possibilità relativamente semplice è la distribuzione di Weibull con P (X> x) = exp [- (x/b)^a].

Adattare i valori come P (X> 1) = 0,1 e P (X> 10) = 0,005, ottengo a = 0,36 eb = 0,1. Ciò implicherebbe che P (X> 40) * 10000 = 1,6, che è un po 'troppo basso, ma P (X> 70) * 10000 = 0,2 che è ragionevole.

EDIT Oh, e per generare una variabile aleatoria Weibull-distribuito da una divisa (0,1) il valore U, basta calcolare b * [- log (1-u)]^(1/a). Questa è la funzione inversa di 1-P (X> x) nel caso in cui ho calcolato male qualcosa.

Answer 2

anni fa, per scritte PHP4, semplicemente scegliere la vostra distribuzione:

<?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

Questo modo ingenuo di farlo probabilmente distorcerà la distribuzione in un modo che non riesco a vedere in questo momento. L'idea è semplicemente di scorrere il tuo primo set di dati, ordinati e come coppie. Quindi randomizzare 15 nuovi numeri tra ciascuna coppia per ottenere il nuovo array.

Esempio rubino, dal momento che non parlo molto PHP. Speriamo che un'idea così semplice dovrebbe essere facile da tradurre in 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

Il modo più semplice (ma non molto efficiente) per generare numeri casuali che seguono una data distribuzione è una tecnica chiamata Von Neumann Rejection.

L'esplosione semplice della tecnica è questa. Crea una scatola che racchiuda completamente la tua distribuzione. (consente di chiamare la distribuzione f) Quindi selezionare un punto casuale (x,y) nella casella. Se y < f(x), quindi utilizzare x come un numero casuale. Se y > f(x), quindi scartare sia x e e scegliere un altro punto. Continua finché non hai una quantità sufficiente di valori da utilizzare. I valori di x che non vengono rifiutati verranno distribuiti in base allo f.