Генератор случайных чисел, который производит степенное распределение? - PullRequest
Новая
       36

Генератор случайных чисел, который производит степенное распределение?

27 голосов
/ 28 мая 2009

Я пишу несколько тестов для приложения Linux для командной строки C ++. Я хотел бы сгенерировать группу целых чисел со степенным распределением / длинным хвостом. То есть я получаю некоторые цифры очень часто, но большинство из них относительно редко.

В идеале это были бы просто магические уравнения, которые я мог бы использовать с rand () или одной из случайных функций stdlib. В противном случае был бы удобен простой в использовании кусок C / C ++.

Спасибо!

Ответы [ 4 ]

34 голосов
/ 28 мая 2009

На этой странице в Wolfram MathWorld обсуждается, как получить степенное распределение из равномерного распределения (именно это обеспечивает большинство генераторов случайных чисел).

Краткий ответ (вывод по вышеуказанной ссылке): 

x = [(x1^(n+1) - x0^(n+1))*y + x0^(n+1)]^(1/(n+1))

, где y - равномерная переменная, n - мощность распределения, x0 и x1 определяют диапазон распределения, и x - ваша распределенная переменная степенного закона.

18 голосов
/ 28 мая 2009

Если вы знаете, какое распределение вы хотите (называемое функцией распределения вероятностей (PDF)) и правильно ли оно нормализовано, вы можете интегрировать его, чтобы получить функцию накопительного распределения (CDF), а затем инвертировать CDF (если это возможно), чтобы получить преобразование, которое вам нужно от равномерного [0,1] распределения до желаемого.

Итак, вы начинаете с определения нужного вам дистрибутива. 

P = F(x)

(для x в [0,1]), затем интегрируется, чтобы дать 

C(y) = \int_0^y F(x) dx

Если это можно перевернуть, вы получите 

y = F^{-1}(C)

Так что назовите rand() и вставьте результат как C в последней строке и используйте y.

Этот результат называется Фундаментальной теоремой выборки. Это хлопотно из-за требования нормализации и необходимости аналитически инвертировать функцию.

В качестве альтернативы вы можете использовать технику отклонения: выбрасывать число равномерно в желаемом диапазоне, затем выбрасывать другое число и сравнивать с PDF в месте, обозначенном вашим первым броском. Отклонить, если второй бросок превышает PDF. Как правило, неэффективно для PDF-файлов с большим количеством областей с низкой вероятностью, например, с длинными хвостами ...

Промежуточный подход включает в себя инвертирование CDF методом грубой силы: вы сохраняете CDF в качестве таблицы поиска и выполняете обратный поиск для получения результата.

Настоящая проблема в том, что простые x^-n распределения не нормируются в диапазоне [0,1], поэтому вы не можете использовать теорему выборки. Попробуйте (x + 1) ^ - n вместо ...

3 голосов
/ 28 октября 2017

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

Одна небольшая разница между формулой Wolfram MathWorld в принятом ответе и другими, возможно, более общими уравнениями заключается в том, что показатель степени степенного закона n (который обычно обозначается как альфа) не имеет явного отрицательного знака. Таким образом, выбранное альфа-значение должно быть отрицательным, обычно от 2 до 3.

x0 и x1 означают нижний и верхний пределы распределения.

Итак, вот оно: 

x1 = 5           # Maximum value
x0 = 0.1         # It can't be zero; otherwise X^0^(neg) is 1/0.
alpha = -2.5     # It has to be negative.
y = runif(1e5)   # Number of samples
x = ((x1^(alpha+1) - x0^(alpha+1))*y + x0^(alpha+1))^(1/(alpha+1))
hist(x, prob = T, breaks=40, ylim=c(0,10), xlim=c(0,1.2), border=F, 
col="yellowgreen", main="Power law density")
lines(density(x), col="chocolate", lwd=1)
lines(density(x, adjust=2), lty="dotted", col="darkblue", lwd=2)

enter image description here

или в логарифмическом масштабе: 

h = hist(x, prob=T, breaks=40, plot=F)
     plot(h$count, log="xy", type='l', lwd=1, lend=2, 
     xlab="", ylab="", main="Density in logarithmic scale")

enter image description here

Вот сводка данных: 

> summary(x)
   Min.   1st Qu.  Median    Mean   3rd Qu.    Max. 
  0.1000  0.1208  0.1584    0.2590  0.2511   4.9388
3 голосов
/ 28 мая 2009

Я не могу комментировать математику, необходимую для получения степенного распределения (другие посты содержат предложения), но я бы посоветовал вам ознакомиться со средствами случайных чисел стандартной библиотеки TR1 C ++ в <random>. Они предоставляют больше функциональных возможностей, чем std::rand и std::srand. Новая система определяет модульный API для генераторов, движков и дистрибутивов и предоставляет набор предустановок.

Включены предустановки распространения:

  • uniform_int
  • bernoulli_distribution
  • geometric_distribution
  • poisson_distribution
  • binomial_distribution
  • uniform_real
  • exponential_distribution
  • normal_distribution
  • gamma_distribution

Когда вы определите свое распределение по степенным законам, вы сможете подключить его к существующим генераторам и двигателям. В книге Пита Бекера Стандартные расширения библиотеки C ++ есть замечательная глава по <random>.

Вот статья о том, как создавать другие дистрибутивы (с примерами для Коши, хи-квадрат, Student t и Snedecor F)

...