Имитация розыгрышей из равномерного распределения с использованием розыгрышей из нормального распределения

Недавно я купил ресурс для интервью с наукой о данных, в котором один из вопросов о вероятности был следующим:

Принимая во внимание ничьи из нормального распределения с известными параметрами, как вы можете имитировать ничьи из равномерного распределения?

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

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

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

Есть мысли о том, как я могу подойти к этому вопросу?

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

1. Обратный метод Бокса-Мюллера : из каждой пары нормалей можно построить две независимые формы: atan2 ( Y , X ) (на интервале [ - π , π ] ) и exp ( - ( X 2 +) Y 2 ) / 2 ) (на интервале [ 0 , 1 ] ).$(X,Y)$$\text{atan2}(Y,X)$$[-\pi, \pi]$$\exp(-(X^2+Y^2)/2)$$[0,1]$

2. Возьмем нормали в группы по два и суммируем их квадраты, чтобы получить последовательность переменных Y 1 , Y 2 , … , Y i , … . Выражения, полученные из пар$\chi^2_2$$Y_1, Y_2, \ldots, Y_i, \ldots$

   $$X_i = \frac{Y_{2i}}{Y_{2i-1}+Y_{2i}}$$

   будет иметь распределение , которое является равномерным.$\text{Beta}(1,1)$

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

3. Поскольку точное распределение коэффициента корреляции Пирсона четырехпарной выборки из стандартного двумерного нормального распределения равномерно распределено по , мы можем просто взять нормали группами по четыре пары (то есть восемь значений в каждый набор) и вернуть коэффициент корреляции этих пар. (Это включает в себя простую арифметику плюс две операции с квадратным корнем.)$[-1,1]$

4. С древних времен было известно, что цилиндрическая проекция сферы (поверхность в трехмерном пространстве) является равной по площади . Это означает, что в проекции равномерного распределения на сфере как горизонтальная координата (соответствующая долготе), так и вертикальная координата (соответствующая широте) будут иметь равномерное распределение. Поскольку стандартное нормальное распределение тривариата сферически симметрично, его проекция на сферу является равномерной. Получение долготы по существу аналогично вычислению угла в методе Бокса-Мюллера ( qv ), но прогнозируемая широта является новой. Проекция на сферу просто нормализует тройку координат и в этой точке z - спроектированная широта. Таким образом, возьмите нормальные переменные в группах по три, X 3 i - 2 , X 3 i - 1 , X 3 i , и вычислите$(x,y,z)$$z$$X_{3i-2}, X_{3i-1}, X_{3i}$

   $$\frac{X_{3i}}{\sqrt{X_{3i-2}^2 + X_{3i-1}^2 + X_{3i}^2}}$$

   для .$i=1, 2, 3, \ldots$

5. Поскольку большинство вычислительных систем представляют собой число в двоичном , равномерное поколение числа , как правило , начинается с получения равномерно распределенные целых чисел между и 2 32 - 1 (или какой - либо высокой мощностью 2 связаны с длиной машинного слова) и перемасштабированием их по мере необходимости. Такие целые числа представлены внутри как строки из 32 двоичных цифр. Мы можем получить независимые случайные биты, сравнивая переменную Normal с ее медианой. Таким образом, достаточно разбить переменные Normal на группы с размером, равным требуемому количеству битов, сравнить каждую из них со средним значением и собрать результирующие последовательности результатов true / false в двоичное число. Написание $0$$2^{32}-1$$2$$32$$k$для числа битов и для знака (то есть H ( x ) = 1, когда x > 0 и H ( x ) = 0 в противном случае) мы можем выразить результирующее нормализованное однородное значение в [ 0 , 1 ) с формулой$H$$H(x)=1$$x\gt 0$$H(x)=0$$[0, 1)$

   $$\sum_{j=0}^{k-1} H(X_{ki - j})2^{-j-1}.$$

   Переменные могут быть взяты из любого непрерывного распределения, медиана которого равна 0 (например, стандартная нормаль); они обрабатываются в группах по k, каждая из которых создает одно такое псевдооднородное значение.$X_n$$0$$k$

6. Отклонение выборки является стандартным, гибким, мощным способом рисования случайных отклонений от произвольных распределений. Предположим, что целевой дистрибутив имеет PDF . Значение Y составлено в соответствии с другим распределением с PDF g . На этапе отклонения равномерное значение U, лежащее между 0 и g ( Y ) , рисуется независимо от Y и сравнивается с f ( Y ) : если оно меньше, $f$$Y$$g$$U$$0$$g(Y)$$Y$$f(Y)$$Y$сохраняется, но в противном случае процесс повторяется. Тем не менее, этот подход кажется круговым: как мы можем сгенерировать однородную переменную с процессом, для начала которого требуется равномерная переменная?

   Ответ заключается в том, что нам на самом деле не нужно унифицировать вариацию для выполнения шага отклонения. Вместо этого (предполагая, что ) мы можем перевернуть справедливую монету, чтобы получить 0 или 1 случайным образом. Это будет интерпретироваться как первый бит в двоичном представлении равномерной переменной U в интервале [ 0 , 1 ) . Когда результат 0 , это означает , что 0 ≤ U < 1 / 2 ; в противном случае, 1 / 2 ≤ U < 1 . $g(Y)\ne 0$$0$$1$$U$$[0,1)$$0$$0 \le U \lt 1/2$$1/2\le U \lt 1$Половина времени, этого достаточно , чтобы решить , шаг: отказ , если , но монета 0 , Y должен быть принят; если F ( Y ) / г ( У ) < 1 / 2 , но монета 1 , Y должна быть отклонена; в противном случае, нам нужно перевернуть монету снова, чтобы получить следующий бит U . Потому что - неважно, какое значение f ( $f(Y)/g(Y) \ge 1/2$$0$$Y$$f(Y)/g(Y) \lt 1/2$$1$$Y$$U$ есть - есть 1 / 2 вероятность остановки после каждого флип, ожидаемое число перестроек только 1 / 2 ( 1 ) + 1 / 4 ( 2 ) + 1 / 8 ( 3 ) + ⋯ + 2 - п ( п ) + ⋯ = 2 .$f(Y)/g(Y)$$1/2$$1/2(1)+1/4(2)+1/8(3)+\cdots+2^{-n}(n)+\cdots=2$

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

   

   Использование Исчисление для оптимизации площади прямоугольника, вы увидите , что ее концы должны лежать на , где его высота равна ехр ( - 1 / +2 ) / $\pm 1$, что делает его площадь чуть больше0,48. Используя эту стандартную нормальную плотность в качествеgиавтоматическиотбрасывая все значения за пределами интервала[-1,1], и в противном случае применяя процедуру отклонения, мы получим равномерные вариации в[-1,1]:$\exp(-1/2)/\sqrt{2\pi}\approx 0.241971$$0.48$$g$$[-1,1]$$[-1,1]$

   • В доле времени нормальное изменение находится за пределами [ - 1 , 1 ] и немедленно отклоняется. ( Φ является стандартным нормальным CDF.)$2\Phi(-1) \approx 0.317$$[-1,1]$$\Phi$

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

   • Общая процедура требует в среднем шагов.$1/(2\exp(-1/2)/\sqrt{2\pi}) \approx 2.07$

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

   $$\sqrt{2 e \pi}\left(1-2\Phi(-1)\right) \approx 2.82137.$$

   Хотя это довольно эффективно, обратите внимание, что (1) вычисление нормального PDF требует вычисления экспоненты и (2) значение должно быть предварительно вычислено раз и навсегда. Это все еще немного меньше, чем метод Бокса-Мюллера ( см. ).$\Phi(-1)$

7. В порядковые статистики равномерного распределения имеют экспоненциальные пробелы. Поскольку сумма квадратов двух нормалей (с нулевым средним) является экспоненциальной, мы можем сгенерировать реализацию независимых форм путем суммирования квадратов пар таких нормалей, вычисления их суммарной суммы и масштабирования результатов, чтобы попасть в интервал. [ 0 , 1 ] и отбрасывая последний (который всегда будет равен 1 ). Это приятный подход, потому что он требует только возведения в квадрат, суммирования и (в конце) одного деления.$n$$[0,1]$$1$

   Значения будут автоматически в порядке возрастания. Если такая сортировка желательна, этот метод в вычислительном отношении превосходит все остальные, поскольку он позволяет избежать затрат O ( n log ( n ) ) на сортировку. Однако, если требуется последовательность независимых униформ, то сортировка этих n значений случайным образом сделает свое дело. Поскольку (как видно из метода Бокса-Мюллера, qv ) отношения каждой пары нормалей не зависят от суммы квадратов каждой пары, у нас уже есть средства для получения этой случайной перестановки: упорядочить кумулятивные суммы по соответствующим отношениям , (Если $n$$O(n\log(n))$$n$$n$очень большой, этот процесс может быть выполнен в меньших группах с небольшой потерей эффективности, так как каждой группе требуется только 2 ( k + 1 ) нормали для создания k однородных значений. Следовательно, для фиксированного k асимптотические вычислительные затраты составляют O ( n log ( k ) ) = O ( n ) , для чего требуется 2 n ( 1 + 1 / k ) нормальных переменных для генерации n однородных значений.)$k$$2(k+1)$$k$$k$$O(n\log(k))$$O(n)$$2n(1+1/k)$$n$

8. В превосходном приближении любой нормальный вариант с большим стандартным отклонением выглядит равномерным в диапазонах с гораздо меньшими значениями. Свернув это распределение в диапазон (взяв только дробные части значений), мы тем самым получим распределение, равномерное для всех практических целей. Это чрезвычайно эффективно, требуя одной из самых простых арифметических операций из всех: просто округлите каждое нормальное значение до ближайшего целого числа и сохраняйте избыток. Простота этого подхода становится убедительной, когда мы исследуем практическую реализацию:$[0,1]$R

   rnorm(n, sd=10) %% 1
   

   надежно создает nоднородные значения в диапазоне за счет просто нормальных изменений и почти без вычислений.$[0,1]$n

   (Даже если стандартное отклонение равно , PDF этого приближения отличается от равномерного PDF, как показано на следующем рисунке, менее чем на одну часть из 10 8 ! Для надежного его обнаружения потребуется образец из 10 16 значений: это уже за пределами возможностей любого стандартного теста на случайность. При большем стандартном отклонении неравномерность настолько мала, что даже не может быть рассчитана. Например, при SD 10, как показано в коде, максимальное отклонение от униформы PDF только 10 - 857. )$1$$10^8$$10^{16}$$10$$10^{-857}$

   

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

Об простоте этих конструкций свидетельствует следующий Rкод, который использует только одну или две строки для большинства из них. Их правильность свидетельствуют в результате ближней равномерной гистограммы на основе независимых величин в каждом случае (требуется около 12 секунд для всех семь моделирования). Для справки - в случае, если вас беспокоит количество вариаций, появляющихся на любом из этих графиков, в конце приводится гистограмма однородных значений, смоделированных с помощью генератора равномерных случайных чисел.$100,000$R

$\chi^2$$1000$$3\%$R

set.seed(17)
n <- 1e5
y <- matrix(rnorm(floor(n/2)*2), nrow=2)
x <- c(atan2(y[2,], y[1,])/(2*pi) + 1/2, exp(-(y[1,]^2+y[2,]^2)/2))
hist(x, main="Box-Mueller")

y <- apply(array(rnorm(4*n), c(2,2,n)), c(3,2), function(z) sum(z^2))
x <- y[,2] / (y[,1]+y[,2])
hist(x, main="Beta")

x <- apply(array(rnorm(8*n), c(4,2,n)), 3, function(y) cor(y[,1], y[,2]))
hist(x, main="Correlation")

n.bits <- 32; x <-  (2^-(1:n.bits)) %*% matrix(rnorm(n*n.bits) > 0, n.bits)
hist(x, main="Binary")

y <- matrix(rnorm(n*3), 3)
x <- y[1, ] / sqrt(apply(y, 2, function(x) sum(x^2)))
hist(x, main="Equal area")

accept <- function(p) { # Using random normals, return TRUE with chance `p`
  p.bit <- x <- 0
  while(p.bit == x) {
    p.bit <- p >= 1/2
    x <- rnorm(1) >= 0
    p <- (2*p) %% 1
  }
  return(x == 0)
}
y <- rnorm(ceiling(n * sqrt(exp(1)*pi/2))) # This aims to produce `n` uniforms
y <- y[abs(y) < 1]
x <- y[sapply(y, function(x) accept(exp((x^2-1)/2)))]
hist(x, main="Rejection")

y <- matrix(rnorm(2*(n+1))^2, 2)
x <- cumsum(y)[seq(2, 2*(n+1), 2)]
x <- x[-(n+1)] / x[n+1]
x <- x[order(y[2,-(n+1)]/y[1,-(n+1)])] 
hist(x, main="Ordered")

x <- rnorm(n) %% 1 # (Use SD of 5 or greater in practice)
hist(x, main="Modular")

x <- runif(n)      # Reference distribution
hist(x, main="Uniform")

$X \sim N(\mu, \sigma^2)$$\Phi_{\mu,\sigma^2}$$\Phi_{\mu,\sigma^2}(X)$$(0,1)$$d \in (0,1)$

$P(\Phi_{\mu,\sigma^2}(X) \leq d) = P(X \leq \Phi_{\mu,\sigma^2}^{-1}(d)) = d$

$d$$1$$d \geq 1$$\Phi_{\mu,\sigma^2}(X)$$(0,1)$$\sigma$$\mathbb{R}$