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

Question

Хорошо, вот моя проблема. Мы смотрим на покупку набора данных от компании, чтобы дополнить наш существующий набор данных. Для целей этого вопроса, скажем, этот набор данных ранжирует места с органическим номером (это означает, что номер, присвоенный одному месту, не имеет отношения к номеру, назначенному другому). Технический диапазон от 0 до бесконечности, но из выборочных наборов, которые я видел, это от 0 до 70. Исходя из выборки, это, безусловно, неравномерное распределение (из 10000, возможно, есть 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, что является разумным.

EDIT Да, и для генерации распределенной по Вейбулю случайной величины из равномерного (0,1) значения U просто вычислите b * [- log (1-u)] ^ (1 / a). Это обратная функция 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

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

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

Answer 4

Этот наивный способ сделать это, скорее всего, искажает дистрибутив так, как я сейчас не вижу. Идея состоит в том, чтобы просто перебрать ваш первый набор данных, отсортированный и в виде пар. Затем рандомизируйте 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!