Генерировать случайные числа с вероятностным распределением

Question

Итак, вот моя проблема. Мы рассматриваем покупку набора данных от компании для увеличения нашего существующего набора данных. Для целей этого вопроса предположим, что этот набор данных занимает места с органическим числом (что означает, что число, присвоенное одному месту, не имеет отношения к числу, присвоенному другому). Технический диапазон от 0 до бесконечности, но из наборов образцов, которые я видел, это от 0 до 70. Исходя из выборки, это, безусловно, не равномерное распределение (из 10 000, возможно, 5 мест со счетом более 40, 50 со счетом более 10 и 1000 со счетом более 1). Прежде чем мы решим купить этот набор, мы хотели бы имитировать его, чтобы мы могли видеть, насколько он может быть полезен.

Итак, чтобы имитировать это, я думал о создании случайного числа для каждого места (около 150 000 случайных чисел). Но я также хочу придерживаться духа данных и держать распределение относительно одинаковым (или, по крайней мере, достаточно близким). Я весь день ломаю голову, пытаясь придумать, как это сделать, и опустел.

Одна мысль, которая у меня была, заключалась в том, чтобы квадрат случайного числа (между 0 и sqrt (70)). Но это будет стоить не менее 1 и более.

Я думаю, что он действительно должен быть гиперболическим в первом квадранте ... Я просто задумываюсь о том, как превратить линейное, равномерное распределение случайных чисел в гиперболическое распределение (если гиперболический - это даже то, что я хочу в первую очередь).

Любые мысли?

Таким образом, чтобы подвести, вот распределение Я хотел бы (приблизительно):

• 40 - 70: 0,02% - 0,05%
• 10 - 40: 0,5% - 1%
• 1 - 10: 10% - 20%
• 0 - 1: Остаток (78,95% - 89,48%)

Answer 1

Посмотрите на распределения, используемые при анализе надежности - они имеют тенденцию иметь эти длинные хвосты. Относительно простая возможность - распределение Вейбулла с P (X> x) = exp [- (x/b)^a].

Применяя ваши значения как P (X> 1) = 0,1 и P (X> 10) = 0,005, я получаю a = 0,36 и b = 0,1. Это означало бы, что P (X> 40) * 10000 = 1,6, что слишком мало, но P (X> 70) * 10000 = 0,2, что является разумным.

РЕДАКТИРОВАТЬ О, и для генерации Вейбуллу-распределенной случайной величины с равномерным (0,1) значением U, просто вычислить B * [- log (1-и)]^(1/а). Это обратная функция 1-P (X> x), если я просчитал что-то.

Answer 2

написано лет назад для PHP4, просто выбрать дистрибутив:

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

Этот наивный способ сделать это, скорее всего, будет искажать распределение в некотором роде, которое я не вижу прямо сейчас. Идея состоит в том, чтобы просто перебрать ваш первый набор данных, отсортированный и как пары. Затем рандомизируйте 15 новых чисел между каждой парой, чтобы получить новый массив.

Пример Ruby, так как я мало говорю о PHP. Надеюсь, такая простая идея должна быть легко переведена на 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

Самый простой (но не очень эффективный) способ генерации случайных чисел, которые следуют за данное распределение является метод, называемый Von Neumann Rejection.

Простое объяснение техники заключается в следующем. Создайте поле, которое полностью охватывает ваш дистрибутив. (давайте позвоним вашему дистрибутиву f) Затем выберите случайную точку (x,y) в поле. Если y < f(x), используйте x как случайное число. Если y > f(x), то отбросьте как x, так и y и выберите другой пункт. Продолжайте, пока не получите достаточное количество значений для использования. Значения x, которые вы не отклоняете, будут распределены в соответствии с f.