Преобразование равномерного распределения в нормальное распределение

106

Как я могу преобразовать равномерное распределение (как и большинство генераторов случайных чисел, например, между 0,0 и 1,0) в нормальное распределение? Что, если я хочу выбрать среднее значение и стандартное отклонение?

Terhorst
источник
3
У вас есть спецификация языка, или это просто общий вопрос алгоритма?
Bill the Lizard
3
Общий вопрос об алгоритме. Мне все равно, на каком языке. Но я бы предпочел, чтобы ответ не полагался на конкретные функции, которые предоставляет только этот язык.
Terhorst

Ответы:

47

Алгоритм Зиккурат довольно эффективен для этого, хотя преобразование Бокса-Мюллера проще реализовать с нуля (а не сумасшедший медленно).

Тайлер
источник
7
Обычные предупреждения о линейных генераторах конгруэнтности относятся к обоим этим методам, поэтому используйте достойный подчиненный генератор. Ура.
dmckee --- котенок экс-модератора
3
Например, Mersenee Twister, или у вас есть другие предложения?
Грегг Линд,
47

Есть много способов:

  • Как не использовать Box Muller. Особенно, если вы рисуете много гауссовских чисел. Box Muller дает результат, который находится в пределах от -6 до 6 (при условии двойной точности. С поплавками ситуация ухудшается). И это действительно менее эффективно, чем другие доступные методы.
  • Ziggurat в порядке, но требует поиска в таблице (и некоторых настроек для конкретной платформы из-за проблем с размером кеша)
  • Отношение форм - мой любимый, только несколько сложений / умножений и логарифм в 1/50 времени (например, посмотрите там ).
  • Обращая CDF является эффективным (и упускать из виду, почему?), То есть быстрые реализации него доступны , если вы поиск Google. Это обязательно для квазислучайных чисел.
Александр К.
источник
2
Вы уверены насчет зажима [-6,6]? Это довольно важный момент, если таковой (и заслуживает упоминания на странице википедии).
redcalx
1
@locster: вот что сказал мне мой учитель (он изучал такие генераторы, и я верю его слову). Возможно, я смогу найти вам ссылку.
Alexandre C.
7
@locster: это нежелательное свойство также характерно для обратного метода CDF. См. Cimat.mx/~src/prope08/randomgauss.pdf . Этого можно избежать, используя равномерный ГСЧ, который имеет ненулевую вероятность дать число с плавающей запятой, очень близкое к нулю. Большинство ГСЧ этого не делают, поскольку они генерируют (обычно 64-битное) целое число, которое затем сопоставляется с [0,1]. Это делает эти методы непригодными для выборки хвостов гауссовских переменных (подумайте о ценообразовании с низким / высоким страйком в вычислительных финансах).
Alexandre C.
6
@AlexandreC. Чтобы прояснить два момента, при использовании 64-битных чисел хвосты выходят либо на 8,57, либо на 9,41 (нижнее значение, соответствующее преобразованию в [0,1) перед взятием журнала). Даже при ограничении до [-6, 6] шансы выйти за пределы этого диапазона составляют примерно 1,98e-9, что достаточно хорошо для большинства людей, даже занимающихся наукой. Для цифр 8,57 и 9,41 это становится 1,04e-17 и 4,97e-21. Эти числа настолько малы, что разница между выборкой Бокса-Мюллера и истинной гауссовой выборкой с точки зрения указанного предела является почти чисто академической. Если вам нужно получше, просто сложите четыре из них и разделите на 2.
CrazyCasta
6
Я думаю, что предложение не использовать преобразование Бокса-Мюллера вводит в заблуждение большой процент пользователей. Приятно знать об ограничении, но, как указывает CrazyCasta, для большинства приложений, которые не сильно зависят от выбросов, вам, вероятно, не нужно об этом беспокоиться. Например, если вы когда-либо полагались на выборку из нормального с помощью numpy, вы полагались на преобразование Бокса-Мюллера (форма полярных координат) github.com/numpy/numpy/blob/… .
Андреас Гривас
30

Изменение распределения одной функции на другую предполагает использование функции, обратной функции, которую вы хотите.

Другими словами, если вы стремитесь к определенной функции вероятности p (x), вы получаете распределение, интегрируя по ней -> d (x) = интеграл (p (x)) и используя его обратное: Inv (d (x)) . Теперь используйте функцию случайной вероятности (которая имеет равномерное распределение) и передайте значение результата через функцию Inv (d (x)). Вы должны получить случайные значения с распределением в соответствии с выбранной вами функцией.

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

Надеюсь, это помогло, и спасибо за небольшое замечание об использовании распределения, а не самой вероятности.

Ади
источник
4
+1 Это очень хорошо работающий метод генерации гауссовских переменных. Обратный CDF может быть эффективно вычислен с помощью метода Ньютона в этом случае (производная равна e ^ {- t ^ 2}), начальное приближение легко получить в виде рациональной дроби, поэтому вам потребуется 3-4 вычисления erf и exp. Это обязательно, если вы используете квазислучайные числа, случай, когда вы должны использовать ровно одно однородное число, чтобы получить гауссово.
Alexandre C.
9
Обратите внимание, что вам нужно инвертировать кумулятивную функцию распределения, а не функцию распределения вероятностей. Александр подразумевает это, но я подумал, что более подробное упоминание об этом не повредит - поскольку ответ, кажется, предполагает PDF-файл
ltjax
Вы можете использовать PDF, если готовы случайным образом выбрать направление относительно среднего; я правильно понимаю?
Марк МакКенна
2
Это называется выборкой с
обратным
1
Вот связанный с этим вопрос в SE с более обобщенным ответом с красивым объяснением.
dashesy 03 окт.15,
23

Вот реализация javascript с использованием полярной формы преобразования Бокса-Мюллера.

/*
 * Returns member of set with a given mean and standard deviation
 * mean: mean
 * standard deviation: std_dev 
 */
function createMemberInNormalDistribution(mean,std_dev){
    return mean + (gaussRandom()*std_dev);
}

/*
 * Returns random number in normal distribution centering on 0.
 * ~95% of numbers returned should fall between -2 and 2
 * ie within two standard deviations
 */
function gaussRandom() {
    var u = 2*Math.random()-1;
    var v = 2*Math.random()-1;
    var r = u*u + v*v;
    /*if outside interval [0,1] start over*/
    if(r == 0 || r >= 1) return gaussRandom();

    var c = Math.sqrt(-2*Math.log(r)/r);
    return u*c;

    /* todo: optimize this algorithm by caching (v*c) 
     * and returning next time gaussRandom() is called.
     * left out for simplicity */
}
user5084
источник
5

Используйте центральную предельную теорему wikipedia entry mathworld entry в своих интересах.

Сгенерируйте n из равномерно распределенных чисел, просуммируйте их, вычтите n * 0,5, и вы получите результат примерно нормального распределения со средним значением, равным 0, и дисперсией, равной (1/12) * (1/sqrt(N))(см. Википедию о равномерном распределении для последнего)

n = 10 дает что-то наполовину приличное быстро. Если вам нужно что-то более чем наполовину приличное, выберите решение Tylers (как указано в статье в Википедии о нормальных дистрибутивах )

Jilles de Wit
источник
1
Это не даст особенно близкого нормального распределения («хвосты» или конечные точки не будут близки к реальному нормальному распределению). Бокс-Мюллер лучше, как предполагали другие.
Питер К.
1
У Бокса Мюллера тоже неправильные хвосты (он возвращает число от -6 до 6 с двойной точностью)
Александр С.
n = 12 (сложить 12 случайных чисел в диапазоне от 0 до 1 и вычесть 6) приводит к stddev = 1 и mean = 0. Затем это можно использовать для генерации любого нормального распределения. Просто умножьте результат на желаемое стандартное отклонение и добавьте среднее значение.
JerryM
3

Я бы использовал Box-Muller. Об этом две вещи:

  1. В итоге вы получаете два значения на итерацию.
    Обычно вы кешируете одно значение и возвращаете другое. При следующем вызове образца вы возвращаете кешированное значение.
  2. Бокс-Мюллер дает Z-оценку.
    Затем вам нужно масштабировать Z-оценку по стандартному отклонению и добавить среднее значение, чтобы получить полное значение в нормальном распределении.
коричнево-коричневый
источник
Как вы масштабируете Z-балл?
Terhorst
3
scaled = mean + stdDev * zScore // дает вам нормальное (среднее значение, stdDev ^ 2)
yoyoyoyosef
2

Где R1, R2 - случайные равномерные числа:

НОРМАЛЬНОЕ РАСПРЕДЕЛЕНИЕ с SD равным 1: sqrt (-2 * log (R1)) * cos (2 * pi * R2)

Это точно ... нет необходимости делать все эти медленные петли!

Эрик Аронести
источник
Прежде чем кто-то поправил меня ... вот приближение, которое я придумал: (1,5- (R1 + R2 + R3)) * 1,88. Мне тоже это нравится.
Эрик Аронести,
2

Кажется невероятным, что я мог бы что-то добавить к этому через восемь лет, но в случае Java я хотел бы указать читателям на метод Random.nextGaussian () , который генерирует для вас гауссово распределение со средним значением 0,0 и стандартным отклонением 1,0.

Простое сложение и / или умножение изменит среднее значение и стандартное отклонение в соответствии с вашими потребностями.

Пепейн Шмитц
источник
1

Стандартный модуль библиотеки Python random имеет то, что вы хотите:

normalvariate (mu, sigma)
Нормальное распределение. mu - среднее значение, а сигма - стандартное отклонение.

Что касается самого алгоритма, взгляните на функцию в random.py в библиотеке Python.

Ручной ввод здесь

Брент.Лонгборо
источник
2
К сожалению, библиотека python использует Kinderman, AJ и Monahan, JF, «Компьютерная генерация случайных величин с использованием отношения равномерных отклонений», ACM Trans Math Software, 3, (1977), pp257-260. Здесь используются две однородные случайные переменные для генерации нормального значения, а не одна, поэтому не очевидно, как использовать ее в качестве сопоставления, которое хотел OP.
Ян
1

Это моя реализация на JavaScript алгоритма P ( полярный метод для нормальных отклонений ) из раздела 3.4.1 книги Дональда Кнута Искусство компьютерного программирования :

function normal_random(mean,stddev)
{
    var V1
    var V2
    var S
    do{
        var U1 = Math.random() // return uniform distributed in [0,1[
        var U2 = Math.random()
        V1 = 2*U1-1
        V2 = 2*U2-1
        S = V1*V1+V2*V2
    }while(S >= 1)
    if(S===0) return 0
    return mean+stddev*(V1*Math.sqrt(-2*Math.log(S)/S))
}
Алессандро Якопсон
источник
0

Я считаю, что вам следует попробовать это в EXCEL : =norminv(rand();0;1). Это произведет произведение случайных чисел, которые должны быть нормально распределены с нулевым средним и объединенной дисперсией. «0» может быть поставлен с любым значением, так что числа будут иметь желаемое среднее, и, изменив «1», вы получите дисперсию, равную квадрату вашего ввода.

Например: =norminv(rand();50;3)будет соответствовать нормально распределенным числам со СРЕДНИМ = 50 РАЗБИРАТЕЛЬНОСТЬ = 9.

Бегемот
источник
0

В Как я могу преобразовать равномерное распределение (как и большинство генераторов случайных чисел, например, между 0,0 и 1,0) в нормальное распределение?

  1. Для программной реализации я знаю пару имен генераторов случайных чисел, которые дают вам псевдооднородную случайную последовательность в [0,1] (Mersenne Twister, Linear Congruate Generator). Назовем это U (x)

  2. Существует математическая область, которая называется теорией вероятностей. Первое: если вы хотите смоделировать rv с интегральным распределением F, вы можете попробовать просто оценить F ^ -1 (U (x)). В пр. Теории было доказано, что такая с.в. будет иметь интегральное распределение F.

  3. Шаг 2 может применяться для генерации rv ~ F без использования каких-либо методов подсчета, когда F ^ -1 может быть получено аналитически без проблем. (например, Exp.distribution)

  4. Чтобы смоделировать нормальное распределение, вы можете вычислить y1 * cos (y2), где y1 ~ равномерно в [0,2pi]. y2 - распределение релей.

В: Что, если мне нужно выбрать среднее и стандартное отклонение по моему выбору?

Вы можете вычислить сигму * N (0,1) + m.

Можно показать, что такое смещение и масштабирование приводят к N (m, сигма)

Брузюз
источник
0

Это реализация Matlab с использованием полярной формы преобразования Бокса-Мюллера :

Функция randn_box_muller.m:

function [values] = randn_box_muller(n, mean, std_dev)
    if nargin == 1
       mean = 0;
       std_dev = 1;
    end

    r = gaussRandomN(n);
    values = r.*std_dev - mean;
end

function [values] = gaussRandomN(n)
    [u, v, r] = gaussRandomNValid(n);

    c = sqrt(-2*log(r)./r);
    values = u.*c;
end

function [u, v, r] = gaussRandomNValid(n)
    r = zeros(n, 1);
    u = zeros(n, 1);
    v = zeros(n, 1);

    filter = r==0 | r>=1;

    % if outside interval [0,1] start over
    while n ~= 0
        u(filter) = 2*rand(n, 1)-1;
        v(filter) = 2*rand(n, 1)-1;
        r(filter) = u(filter).*u(filter) + v(filter).*v(filter);

        filter = r==0 | r>=1;
        n = size(r(filter),1);
    end
end

И ссылаясь на histfit(randn_box_muller(10000000),100);это результат: Box-Muller Matlab Histfit

Очевидно, это действительно неэффективно по сравнению со встроенным в Matlab randn .

Madx
источник
0

У меня есть следующий код, который может помочь:

set.seed(123)
n <- 1000
u <- runif(n) #creates U
x <- -log(u)
y <- runif(n, max=u*sqrt((2*exp(1))/pi)) #create Y
z <- ifelse (y < dnorm(x)/2, -x, NA)
z <- ifelse ((y > dnorm(x)/2) & (y < dnorm(x)), x, z)
z <- z[!is.na(z)]
great_minds_think_alike
источник
0

Также проще использовать реализованную функцию rnorm (), поскольку она быстрее, чем написание генератора случайных чисел для нормального распределения. См. Следующий код как доказательство

n <- length(z)
t0 <- Sys.time()
z <- rnorm(n)
t1 <- Sys.time()
t1-t0
Питервитетбитер
источник
-2
function distRandom(){
  do{
    x=random(DISTRIBUTION_DOMAIN);
  }while(random(DISTRIBUTION_RANGE)>=distributionFunction(x));
  return x;
}

источник
Однако не гарантировано возвращение, не так ли? ;-)
Питер К.
5
Случайные числа слишком важны, чтобы их оставлять на волю случая.
Дрю Ноукс
Не отвечает на вопрос - нормальное распределение имеет бесконечную область.
Мэтт