Генерация случайных чисел в соответствии с нормальным распределением в C / C ++

114

Как я могу легко генерировать случайные числа в соответствии с нормальным распределением в C или C ++?

Я не хочу использовать Boost.

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

Дэмиен
источник
2
Дубликат одного или другого из stackoverflow.com/questions/75677/… и stackoverflow.com/questions/1109446/…
dmckee --- котенок бывшего модератора

Ответы:

92

Существует множество методов генерации чисел с распределением по Гауссу из обычного ГСЧ .

Преобразование Бокса-Мюллера обычно используется. Он правильно производит значения с нормальным распределением. Математика проста. Вы генерируете два (однородных) случайных числа, и применяя к ним формулу, вы получаете два нормально распределенных случайных числа. Верните один, а другой сохраните для следующего запроса случайного числа.

С. Лотт
источник
10
Если вам нужна скорость, то полярный метод быстрее. А алгоритм Зикгурата даже больше (правда, писать намного сложнее).
Джоуи,
2
нашел здесь реализацию зиккурата people.sc.fsu.edu/~jburkardt/c_src/ziggurat/ziggurat.html Это довольно полно.
dwbrito
24
Обратите внимание, что в C ++ 11 добавлено то, std::normal_distributionчто делает именно то, что вы просите, не вдаваясь в математические детали.
3
std :: normal_distribution не гарантирует единообразия на всех платформах. Сейчас я провожу тесты, и MSVC предоставляет другой набор значений, чем, например, Clang. Механизмы C ++ 11, похоже, генерируют одни и те же последовательности (при одном и том же начальном значении), но дистрибутивы C ++ 11, похоже, реализованы с использованием разных алгоритмов на разных платформах.
Arno Duvenhage
47

C ++ 11

С ++ 11 предлагает std::normal_distribution, и я бы пошел именно по этому пути.

C или более ранний C ++

Вот несколько решений в порядке возрастания сложности:

  1. Добавьте 12 одинаковых случайных чисел от 0 до 1 и вычтите 6. Это будет соответствовать среднему значению и стандартному отклонению нормальной переменной. Очевидным недостатком является то, что диапазон ограничен до ± 6 - в отличие от истинного нормального распределения.

  2. Преобразование Бокса-Мюллера. Это перечислено выше и относительно просто реализовать. Однако, если вам нужны очень точные образцы, имейте в виду, что преобразование Бокса-Мюллера в сочетании с некоторыми однородными генераторами страдает аномалией, называемой эффектом Neave 1 .

  3. Для большей точности я предлагаю нарисовать униформу и применить обратное кумулятивное нормальное распределение, чтобы получить нормально распределенные переменные. Вот очень хороший алгоритм для обратных кумулятивных нормальных распределений.

1. HR Neave, "Об использовании преобразования Бокса-Мюллера с мультипликативными конгруэнтными генераторами псевдослучайных чисел", Applied Statistics, 22, 92-97, 1973

Питер Г.
источник
Может быть, у вас есть еще одна ссылка на pdf-файл об эффекте Нива? или ссылка на исходную статью журнала? спасибо
pyCthon
2
@stonybrooknick Добавлена ​​исходная ссылка. Прикольное замечание: во время поиска в Google "box muller neave", чтобы найти ссылку, этот самый вопрос о stackoverflow возник на первой странице результатов!
Питер Г.
да, это не все хорошо известны за пределами определенных небольших сообществ и групп по интересам
pyCthon
@Peter G. Почему кто-то проголосовал против вашего ответа? - возможно, тот же человек оставил и мой комментарий ниже, что меня устраивает, но я думаю, что ваш ответ был очень хорошим. Было бы хорошо, если бы так, чтобы голоса против приводили к реальным комментариям ... Я подозреваю, что большинство голосов против старых тем просто легкомысленно и трогательно.
Pete855217
«Сложите 12 одинаковых чисел от 0 до 1 и вычтите 6.» - распределение этой переменной будет иметь нормальное распределение? Можете ли вы предоставить ссылку на вывод, потому что во время вывода центральной предельной теоремы n -> + inf является очень необходимым предположением.
bruziuz
31

Быстрый и простой способ - просто просуммировать несколько равномерно распределенных случайных чисел и взять их среднее значение. См. Центральную предельную теорему для полного объяснения того, почему это работает.

Пол Р
источник
+1 Очень интересный подход. Проверено, действительно ли дает нормально распределенные подгруппы для небольших групп?
Morlock
4
@Morlock Чем больше вы усредняете выборок, тем ближе вы подходите к распределению Гаусса. Если ваше приложение предъявляет строгие требования к точности распределения, вам может быть лучше использовать что-то более строгое, например Box-Muller, но для многих приложений, например генерации белого шума для аудиоприложений, вы можете обойтись довольно небольшим числом усредненных выборок (например, 16).
Paul R
2
Кроме того, как вы параметризуете это, чтобы получить определенную дисперсию, скажем, вам нужно среднее значение 10 со стандартным отклонением 1?
Morlock,
1
@Ben: не могли бы вы указать мне на эффективный алгоритм для этого? Я использовал только технику усреднения для генерации приблизительно гауссовского шума для обработки звука и изображения с ограничениями в реальном времени - если есть способ добиться этого за меньшее количество тактовых циклов, это может быть очень полезно.
Paul R
1
@Petter: вы, наверное, правы в общем случае для значений с плавающей запятой. Тем не менее, есть области приложений, такие как аудио, где вам нужен быстрый гауссовский шум с целым числом (или с фиксированной точкой), и точность не слишком важна, где простой метод усреднения более эффективен и полезен (особенно для встроенных приложений, где может даже не быть быть аппаратной поддержкой с плавающей запятой).
Paul R
24

Я создал проект с открытым исходным кодом C ++ для теста генерации нормально распределенных случайных чисел .

Он сравнивает несколько алгоритмов, в том числе

  • Метод центральной предельной теоремы
  • Преобразование Бокса-Мюллера
  • Полярный метод Марсальи
  • Алгоритм зиккурата
  • Метод обратного преобразования выборки.
  • cpp11randomиспользует C ++ 11 std::normal_distributionwith std::minstd_rand(на самом деле это преобразование Бокса-Мюллера в clang).

Результаты floatверсии single-precision ( ) на iMac Corei5-3330S@2.70GHz, clang 6.1, 64-bit:

normaldistf

Для правильности программа проверяет среднее значение, стандартное отклонение, асимметрию и эксцесс образцов. Было обнаружено, что метод CLT путем суммирования 4, 8 или 16 однородных чисел не имеет хорошего эксцесса, как другие методы.

Алгоритм зиккурата имеет лучшую производительность, чем другие. Однако он не подходит для параллелизма SIMD, так как требует поиска в таблице и переходов. Box-Muller с набором инструкций SSE2 / AVX намного быстрее (x1,79, x2,99), чем не-SIMD версия алгоритма зиккурата.

Поэтому я предлагаю использовать Box-Muller для архитектуры с наборами инструкций SIMD, в противном случае это может быть зиккуратным.


PS тест использует простейший LCG PRNG для генерации равномерно распределенных случайных чисел. Поэтому для некоторых приложений этого может быть недостаточно. Но сравнение производительности должно быть справедливым, потому что все реализации используют один и тот же PRNG, поэтому тест в основном проверяет производительность преобразования.

Майло Ип
источник
2
«Но сравнение производительности должно быть справедливым, потому что все реализации используют один и тот же ГПСЧ» .. За исключением того, что BM использует один входной RN на выход, тогда как CLT использует намного больше и т.д ... поэтому время для генерации однородного случайного # имеет значение.
greggo
14

Вот пример C ++, основанный на некоторых ссылках. Это быстро и грязно, лучше не изобретать заново и не использовать библиотеку ускорения.

#include "math.h" // for RAND, and rand
double sampleNormal() {
    double u = ((double) rand() / (RAND_MAX)) * 2 - 1;
    double v = ((double) rand() / (RAND_MAX)) * 2 - 1;
    double r = u * u + v * v;
    if (r == 0 || r > 1) return sampleNormal();
    double c = sqrt(-2 * log(r) / r);
    return u * c;
}

Вы можете использовать график QQ, чтобы изучить результаты и посмотреть, насколько хорошо оно приближается к реальному нормальному распределению (ранжируйте ваши образцы 1..x, превратите ранги в пропорции от общего количества x, т.е. сколько образцов, получите z-значения и начертите их. Желаемый результат - прямая линия вверх).

Пит855217
источник
1
Что такое sampleNormalManual ()?
Решение Пазлов
@solvingPuzzles - извините, поправил код. Это рекурсивный вызов.
Pete855217
1
Это обязательно приведет к сбою в каком-то редком событии (демонстрация приложения вашему боссу - это звоночек?). Это должно быть реализовано с помощью цикла, а не с использованием рекурсии. Метод выглядит незнакомым. Что такое источник / как он называется?
the swine
Бокс-Мюллер переписан из java-реализации. Как я уже сказал, это быстро и грязно, не стесняйтесь исправлять.
Pete855217
1
FWIW, многие компиляторы смогут превратить этот конкретный рекурсивный вызов в «переход к началу функции». Вопрос в том, хотите ли вы на это рассчитывать :-) Кроме того, вероятность того, что потребуется> 10 итераций, равна 1 из 4,8 миллиона. p (> 20) - квадрат этого и т. д.
greggo
12

Используйте std::tr1::normal_distribution.

Пространство имен std :: tr1 не является частью boost. Это пространство имен, которое содержит дополнения к библиотеке из Технического отчета C ++ 1 и доступно в современных компиляторах Microsoft и gcc, независимо от boost.

JoeG
источник
25
Он не просил стандарта, он просил «не повышать».
JoeG
12

Вот как вы создаете образцы на современном компиляторе C ++.

#include <random>
...
std::mt19937 generator;
double mean = 0.0;
double stddev  = 1.0;
std::normal_distribution<double> normal(mean, stddev);
cerr << "Normal: " << normal(generator) << endl;
Петтер
источник
generatorдействительно должны быть посеяны.
Уолтер
Он всегда засевается. Есть семя по умолчанию.
Петтер
4

Если вы используете C ++ 11, вы можете использовать std::normal_distribution:

#include <random>

std::default_random_engine generator;
std::normal_distribution<double> distribution(/*mean=*/0.0, /*stddev=*/1.0);

double randomNumber = distribution(generator);

Есть много других распределений, которые можно использовать для преобразования вывода механизма случайных чисел.

Дрю Ноукс
источник
Об этом уже упоминал Бен ( stackoverflow.com/a/11977979/635608 )
Mat
3

Я следил за определением PDF-файла, приведенным в http://www.mathworks.com/help/stats/normal-distribution.html, и пришел к следующему:

const double DBL_EPS_COMP = 1 - DBL_EPSILON; // DBL_EPSILON is defined in <limits.h>.
inline double RandU() {
    return DBL_EPSILON + ((double) rand()/RAND_MAX);
}
inline double RandN2(double mu, double sigma) {
    return mu + (rand()%2 ? -1.0 : 1.0)*sigma*pow(-log(DBL_EPS_COMP*RandU()), 0.5);
}
inline double RandN() {
    return RandN2(0, 1.0);
}

Возможно, это не лучший подход, но он довольно простой.

MJVC
источник
-1 Не работает, например, для RANDN2 (0.0, d + 1.0). Макросы печально известны этим.
Петтер
Макрос будет ошибкой , если rand()из RANDUвозвращает ноль, так как Ln (0) не определено.
InterDist
Вы действительно пробовали этот код? Похоже, вы создали функцию, которая генерирует числа с распределением Рэлея . Сравните с преобразованием Бокса – Мюллера , где они умножаются на cos(2*pi*rand/RAND_MAX), а вы умножаете на (rand()%2 ? -1.0 : 1.0).
HelloGoodbye
1

Список часто задаваемых вопросов на comp.lang.c содержит три различных способа простого генерирования случайных чисел с распределением Гаусса.

Вы можете взглянуть на это: http://c-faq.com/lib/gaussian.html

Delgan
источник
1

Реализация Бокса-Мюллера:

#include <cstdlib>
#include <cmath>
#include <ctime>
#include <iostream>
using namespace std;
 // return a uniformly distributed random number
double RandomGenerator()
{
  return ( (double)(rand()) + 1. )/( (double)(RAND_MAX) + 1. );
}
 // return a normally distributed random number
double normalRandom()
{
  double y1=RandomGenerator();
  double y2=RandomGenerator();
  return cos(2*3.14*y2)*sqrt(-2.*log(y1));
}

int main(){
double sigma = 82.;
double Mi = 40.;
  for(int i=0;i<100;i++){
double x = normalRandom()*sigma+Mi;
    cout << " x = " << x << endl;
  }
  return 0;
}
Сисадмин
источник
1

Существуют различные алгоритмы обратного кумулятивного нормального распределения. Самые популярные в количественном финансировании тестируются на http://chasethedevil.github.io/post/monte-carlo--inverse-cumulative-normal-distribution/

На мой взгляд, существует не так много стимулов для использования что - то другое , чем алгоритм AS241 от Wichura : это машина точность, надежность и быстро. Узкие места редко возникают при генерации гауссовских случайных чисел.

Кроме того, он показывает недостатки зиккуратоподобных подходов.

Главный ответ здесь выступает за Box-Müller, вы должны знать, что у него есть известные недостатки. Цитирую https://www.sciencedirect.com/science/article/pii/S0895717710005935 :

в литературе Бокса – Мюллера иногда считают немного хуже, в основном по двум причинам. Во-первых, если применить метод Бокса – Мюллера к числам из плохого линейного конгруэнтного генератора, преобразованные числа обеспечат очень плохое покрытие пространства. Графики преобразованных чисел со спиралевидными хвостами можно найти во многих книгах, в первую очередь в классической книге Рипли, который, вероятно, первым сделал это наблюдение "

Джерек
источник
0

1) Графически интуитивно понятный способ генерации гауссовских случайных чисел - это использование чего-то похожего на метод Монте-Карло. Вы должны сгенерировать случайную точку в рамке вокруг кривой Гаусса, используя генератор псевдослучайных чисел в C. Вы можете вычислить, находится ли эта точка внутри или под распределением Гаусса, используя уравнение распределения. Если эта точка находится внутри распределения Гаусса, то у вас есть гауссовское случайное число в качестве значения x точки.

Этот метод не идеален, потому что технически кривая Гаусса стремится к бесконечности, и вы не можете создать коробку, которая приближается к бесконечности в измерении x. Но кривая Гуасса довольно быстро приближается к 0 в измерении y, так что я не буду об этом беспокоиться. Ограничение размера ваших переменных в C может быть более ограничивающим фактором для вашей точности.

2) Другой способ - использовать Центральную предельную теорему, которая утверждает, что при добавлении независимых случайных величин они образуют нормальное распределение. Помня об этой теореме, вы можете аппроксимировать гауссовское случайное число, добавив большое количество независимых случайных величин.

Эти методы не самые практичные, но этого следует ожидать, если вы не хотите использовать уже существующую библиотеку. Имейте в виду, что этот ответ исходит от человека, практически не имеющего опыта в расчетах или статистике.

дан дан
источник
0

Метод Монте-Карло . Наиболее интуитивно понятный способ сделать это - использовать метод Монте-Карло. Возьмите подходящий диапазон -X, + X. Большие значения X приведут к более точному нормальному распределению, но для схождения потребуется больше времени. а. Выберите случайное число z от -X до X. b. Сохраняйте с вероятностью, N(z, mean, variance)где N - гауссово распределение. В противном случае отбросьте и вернитесь к шагу (а).

Джагат
источник
-3

Компьютер - детерминированное устройство. В расчетах нет случайности. Более того, арифметическое устройство в CPU может вычислять сумму по некоторому конечному набору целых чисел (выполняя вычисление в конечном поле) и конечному набору действительных рациональных чисел. А также выполнял побитовые операции. Математика имеет дело с более хорошими наборами, такими как [0.0, 1.0] с бесконечным количеством баллов.

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

Существуют алгоритмы, которые называются - генератор псевдослучайных чисел. Как я понял, цель псевдослучайного генератора - имитировать случайность. И критерий добродетели: - эмпирическое распределение сходится (в некотором смысле - поточечное, равномерное, L2) к теоретическому - значения, которые вы получаете от случайного генератора, кажутся независимыми. Конечно, это неправда с «реальной точки зрения», но мы предполагаем, что это правда.

Один из популярных методов - вы можете суммировать 12 irv с равномерными распределениями .... Но, честно говоря, при выводе центральной предельной теоремы с помощью преобразования Фурье, рядов Тейлора необходимо пару раз иметь n -> + inf предположений. Так, например, теоретически - лично я не понимаю, как люди выполняют суммирование 12 irv с равномерным распределением.

В университете я изучал теорию вероятностей. И, в частности, для меня это просто математический вопрос. В университете я увидел такую ​​модель:


double generateUniform(double a, double b)
{
  return uniformGen.generateReal(a, b);
}

double generateRelei(double sigma)
{
  return sigma * sqrt(-2 * log(1.0 - uniformGen.generateReal(0.0, 1.0 -kEps)));
}
double generateNorm(double m, double sigma)
{
  double y2 = generateUniform(0.0, 2 * kPi);
  double y1 = generateRelei(1.0);
  double x1 = y1 * cos(y2);
  return sigma*x1 + m;
}

Такой способ, как это сделать, был всего лишь примером, я думаю, существуют другие способы его реализации.

Доказательство ее правильности можно найти в книге "Москва, БМГТУ, 2004: XVI Теория вероятностей, пример 6.12, с.246-247" в издании Крищенко Александра Петровича ISBN 5-7038-2485-0

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

Брузюз
источник
У меня несколько голосов против. Дайте мне знать, что здесь плохого?
bruziuz
Вопрос в том, как сгенерировать псевдослучайные числа на компьютере (я знаю, здесь нет слов), это не вопрос математического существования.
user2820579
Да, ты прав. И ответ заключается в том, как сгенерировать псевдослучайное число с нормальным распределением на основе генератора, который имеет равномерное распределение. Исходный код предоставлен, вы можете переписать его на любом языке.
bruziuz
Конечно, я думаю, что этот парень ищет, например, «Числовые рецепты на C / C ++». Между прочим, чтобы дополнить наше обсуждение, авторы этой последней книги дают интересные ссылки на пару псевдослучайных генераторов, которые соответствуют стандартам «приличных» генераторов.
user2820579
1
Я сделал резервную копию здесь: sites.google.com/site/burlachenkok/download
bruziuz