Генератор псевдослучайных чисел - экспоненциальное распределение

Я хотел бы создать несколько псевдослучайных чисел, и до сих пор я был очень доволен библиотекой.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. Пусть будет 1/rate. Установите на 0.
  2. Сгенерируйте однородную (0,) случайную переменную (например, NextDouble()*scale), назови это .
  3. Задавать val на и установите на 1.
  4. Сгенерируйте однородную (0,) случайную переменную, назовите ее.
  5. Если больше, установите u к v, затем установите 1 минус accept, затем переходите к шагу 4.
  6. Если accept равен 1, вернуть val + highpart.
  7. Добавлять 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))
}

Конечно, это формула возведения в квадрат случайного числа, поэтому вы генерируете случайное число вдоль квадратичной кривой.

Другие вопросы по тегам