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

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

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

Обходной путь для suhosin.mt_srand.ignore для последовательного перетасовки массива в PHP?

Я думал, что это было квадратное случайное число (между 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%)

Посмотрите на распределения, используемые в анализе надежности – они имеют тенденцию иметь эти длинные хвосты. Относительно простая возможность – распределение Вейбулла с 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 Oh, и для генерации случайной величины, распределенной по Вейбулу, из равномерного (0,1) значения U, просто вычислим b * [- log (1-u)] ^ (1 / a). Это обратная функция 1-P (X> x), если я просчитал что-то.

Письменные годы назад для 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; } ?>

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

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

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