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

Question

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

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

Спасибо!

Answer 1

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

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

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

где й является однородной вариацией, н это распределительная мощность, х0 и х1 определить диапазон распределения и Икс это ваша степенная распределенная переменная.

Answer 2

Если вы знаете желаемое распределение (называемое функцией распределения вероятностей (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...

Answer 3

Я не могу комментировать математические вычисления, необходимые для получения степенного распределения (в других сообщениях есть предложения), но я бы посоветовал вам ознакомиться со средствами случайных чисел стандартной библиотеки 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>.

Вот статья о том, как создавать другие распределения (с примерами для Коши, Хи-квадрат, Стьюдента t и Snedecor F)

Answer 4

Я просто хотел провести реальное моделирование в дополнение к (справедливо) принятому ответу.Хотя в 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)

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

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")

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

> summary(x)
   Min.   1st Qu.  Median    Mean   3rd Qu.    Max. 
  0.1000  0.1208  0.1584    0.2590  0.2511   4.9388