Генератор псевдослучайных чисел - экспоненциальное распределение
Я хотел бы создать несколько псевдослучайных чисел, и до сих пор я был очень доволен библиотекой.Net Random.Next(int min, int max)
функция. Предполагается, что PRNG этого сорта используют равномерное распределение, но я бы очень хотел сгенерировать некоторые числа с использованием экспоненциального распределения.
Я программирую на C#, хотя я приму псевдокод или C++, Java или тому подобное.
Любые предложения / фрагменты кода / алгоритмы / мысли?
9 ответов
Поскольку у вас есть доступ к унифицированному генератору случайных чисел, генерировать случайное число, распределенное по другому распределению, CDF которого вы знаете, легко, используя метод инверсии.
Итак, сгенерируйте равномерное случайное число, u
, в [0,1)
, а затем рассчитать x
от:
x = log(1-u)/(
-λ)
,
где λ - показатель скорости экспоненциального распределения. Сейчас, x
случайное число с экспоненциальным распределением. Обратите внимание, что log
выше ln
натуральный логарифм.
Фундаментальная теорема выборки гласит, что если вы можете нормализовать, интегрировать и инвертировать желаемое распределение, вы свободны дома.
Если у вас есть желаемое распределение F(x)
нормализовано на [a,b]
, Вы вычисляете
C(y) = \int_a^y F(x) dx
инвертировать это, чтобы получить C^{-1}
бросить z
равномерно на [0,1) и находим
x_i = C^{-1}(z_i)
который будет иметь желаемое распределение.
В твоем случае: F(x) = ke^{-kx}
и я буду считать, что вы хотите [0,infinity]
, Мы получаем:
C(y) = 1 - e^{-ky}
который обратим, чтобы дать
x = -1/k ln(1 - z)
для Z брошен равномерно на [0,1)
,
Но, честно говоря, использование хорошо отлаженной библиотеки разумнее, если вы не делаете это для своего собственного назидания.
Это формула, которую я нашел в Википедии:
Т = -Ln (и) / λ
Мы создаем случайное число с равномерным распределением (u) в [0,1] и получаем x:
Случайный R = новый Случайный ();
двойной u = R. NextDouble();
двойной x = -Math.Log (u) / (λ);
Если вам нужны хорошие случайные числа, рассмотрите возможность ссылки на подпрограммы gsl: http://www.gnu.org/software/gsl/. У них есть рутина gsl_ran_exponential
, Если вы хотите генерировать случайные числа, используя встроенный генератор с равномерным распределением на [0, 1) (например, u=Random.Next(0, N-1)/N, для некоторого большого N), тогда просто используйте:
-mu * log (1-u)
Смотрите randist/exponential.c в источнике gsl.
РЕДАКТИРОВАТЬ: просто для сравнения с некоторыми более поздними ответами - это эквивалентно с mu = 1/ лямбда. mu - это среднее значение распределения, также называемое параметром масштаба на странице википедии, с которым связан OP, а лямбда - параметр скорости.
Одно интересное свойство экспоненциального распределения: рассмотрим процесс прибытия с экспоненциальным временем взаимодействия. Возьмите любой период времени (t1, t2) и прибытия в этот период. Эти прибытия равномерно распределены между t1 и t2. Шелдон Росс. Стохастические процессы.
Если у меня есть генератор псевдослучайных чисел и по какой-то причине (например, мое программное обеспечение не может вычислять журналы), вы не хотите выполнять вышеуказанное преобразование, но хотите получить экспоненциальное значение rv со средним значением 1,0.
Вы можете:
1) Создайте 1001 U(0,1) случайных величин.
2) Сортировать по порядку
3) Вычтите второе из первого, третье из второго, чтобы получить 1000 отличий.
4) Эти различия являются экспоненциальными RV с распределением со средним значением = 1,0.
Думаю, менее эффективным, но средством для достижения той же цели.
Библиотека с открытым исходным кодом Uncommons Maths от Dan Dyer предоставляет генераторы случайных чисел, вероятностные распределения, комбинаторику и статистику для Java.
Среди других ценных классов, ExponentialGenerator
по существу реализовал идею, объясненную @Alok Singhal. В его учебном блоге приведен фрагмент кода для имитации случайного события, которое происходит в среднем 10 раз в минуту:
final long oneMinute = 60000;
Random rng = new MersenneTwisterRNG();
// Generate events at an average rate of 10 per minute.
ExponentialGenerator gen = new ExponentialGenerator(10, rng);
boolean running = true;
while (true)
{
long interval = Math.round(gen.nextValue() * oneMinute);
Thread.sleep(interval);
// Fire event here.
}
Конечно, если вы предпочитаете единицу времени per second
(вместо a minute
здесь), вам просто нужно установить final long oneMinute = 1000
,
Углубляясь в исходный код метода nextValue()
из ExponentialGenerator
вы найдете так называемую выборку обратного преобразования, описанную в Generating_exponential_variates [wiki]:
public Double nextValue()
{
double u;
do
{
// Get a uniformly-distributed random double between
// zero (inclusive) and 1 (exclusive)
u = rng.nextDouble();
} while (u == 0d); // Reject zero, u must be positive for this to work.
return (-Math.log(u)) / rate.nextValue();
}
PS: недавно я использую библиотеку Uncommons Maths. Спасибо Дэн Дайер.
Есть еще один способ сгенерировать экспоненту (
rate
) случайная величина, хотя в настоящее время это не так удобно, как использование логарифмов. Он основан на алгоритме Джона фон Неймана (1951) и использует только сравнения.
- Пусть будет
1/rate
. Установите на 0. - Сгенерируйте однородную (0,) случайную переменную (например,
NextDouble()*scale
), назови это . - Задавать
val
на и установите на 1. - Сгенерируйте однородную (0,) случайную переменную, назовите ее.
- Если больше, установите
u
кv
, затем установите 1 минусaccept
, затем переходите к шагу 4. - Если accept равен 1, вернуть
val + highpart
. - Добавлять
scale
кhighpart
и переходите к шагу 2.
ИСПОЛЬЗОВАННАЯ ЛИТЕРАТУРА:
- Фон Нейман, Дж., «Различные методы, используемые в связи со случайными цифрами», 1951.
Если я понимаю вашу проблему, и вы можете принять конечное число PRNG, вы можете использовать такой подход, как:
- Создайте массив, в котором каждый элемент находится в вашем экспоненциальном распределении
- Создайте PRNG, который является целочисленным индексом в массиве. Вернуть элемент в массиве по этому индексу.
Это то, что я использовал, когда столкнулся с похожими требованиями:
// sorry.. pseudocode, mine was in Tcl:
int weighted_random (int max) {
float random_number = rand();
return floor(max - ceil( max * random_number * random_number))
}
Конечно, это формула возведения в квадрат случайного числа, поэтому вы генерируете случайное число вдоль квадратичной кривой.