Generar números aleatorios con distribución probabilística

Question

Ok, este es mi problema. Estamos buscando comprar un conjunto de datos de una empresa para aumentar nuestro conjunto de datos existente. A los efectos de esta pregunta, digamos que este conjunto de datos clasifica lugares con un número orgánico (lo que significa que el número asignado a un lugar no tiene relación con el número asignado a otro). El rango técnico es de 0 hasta el infinito, pero de los conjuntos de muestras que he visto, es de 0 a 70. Según la muestra, definitivamente no es una distribución uniforme (de 10,000 hay tal vez 5 lugares con una puntuación superior a 40, 50 con un puntaje mayor a 10 y 1000 con un puntaje mayor a 1). Antes de decidir comprar este conjunto, nos gustaría simularlo para que podamos ver lo útil que puede ser.

Por lo tanto, para simularlo, he estado pensando en generar un número aleatorio para cada lugar (aproximadamente 150,000 números aleatorios). Pero también quiero mantener el espíritu de los datos y mantener la distribución relativamente igual (o al menos razonablemente cercana). He estado trabajando en mi cerebro todo el día tratando de pensar en una forma de hacerlo, y he salido vacío.

Un pensamiento que tuve fue cuadrar el número aleatorio (entre 0 y sqrt (70)). Pero eso favorecería tanto a menos de 1 como a números más grandes.

Estoy pensando que la distribución real debería ser hiperbólica en el primer cuadrante ... Solo me estoy quedando sin información sobre cómo convertir una distribución lineal y uniforme de números aleatorios en una distribución hiperbólica (si es incluso hiperbólica) querer en primer lugar).

¿Alguna idea?

tanto, para resumir, aquí está la distribución me gustaría (aproximadamente):

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

Answer 1

Observe las distribuciones utilizadas en el análisis de confiabilidad: tienden a tener estas colas largas. Una posibilidad relativamente simple es la distribución de Weibull con P (X> x) = exp [- (x/b)^a].

Ajustando sus valores como P (X> 1) = 0.1 y P (X> 10) = 0.005, obtengo a = 0.36 yb = 0.1. Esto implicaría que P (X> 40) * 10000 = 1.6, que es un poco demasiado bajo, pero P (X> 70) * 10000 = 0.2 que es razonable.

EDITAR Oh, y para generar una variable aleatoria distribuida-Weibull de un uniforme (0,1) valor U, simplemente calcular b * [- log (1-u)]^(1/a). Esta es la función inversa de 1-P (X> x) en caso de que haya calculado mal algo.

Answer 2

hace

años escrito por PHP4, simplemente escoger su distribución:

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

Esta forma ingenua de hacerlo muy probablemente sesgará la distribución de alguna manera que no puedo ver en este momento. La idea es simplemente iterar sobre su primer conjunto de datos, ordenados y como pares. Luego aleatorice 15 números nuevos entre cada par para obtener la nueva matriz.

Ejemplo de Ruby, ya que no hablo mucho PHP. Con suerte, una idea tan simple debería ser fácil de traducir a 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

La forma más sencilla (pero no muy eficiente) para generar números aleatorios que siguen una distribución dada es una técnica llamada Von Neumann Rejection.

La explicación simple de la técnica es esta. Crea una caja que encierre completamente tu distribución. (vamos a llamar a su distribución f) A continuación, elija un punto aleatorio (x,y) en el cuadro. Si es y < f(x), utilice x como número aleatorio. Si es y > f(x), deseche ambos x y y y elija otro punto. Continúe hasta que tenga una cantidad suficiente de valores para usar. Los valores de x que no rechace se distribuirán según f.