Как генерировать коррелированные случайные числа (с учетом средних, дисперсий и степени корреляции)?

Извините, если это кажется слишком основополагающим, но я думаю, что я просто пытаюсь подтвердить понимание здесь. У меня есть чувство, что я должен сделать это в два этапа, и я начал пытаться получить матрицы корреляции, но это только начинает казаться действительно вовлеченным. Я ищу краткое объяснение (в идеале с подсказками для решения псевдокода) хорошего, идеально быстрого способа генерирования коррелированных случайных чисел.

Учитывая две псевдослучайные переменные роста и веса с известными средними и дисперсиями, а также данную корреляцию, я думаю, что я в основном пытаюсь понять, как должен выглядеть этот второй шаг:

   height = gaussianPdf(height.mean, height.variance)
   weight = gaussianPdf(correlated_mean(height.mean, correlation_coefficient), 
                        correlated_variance(height.variance, 
                        correlation_coefficient))

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

Чтобы ответить на ваш вопрос о «хорошем, в идеале быстром способе генерирования коррелированных случайных чисел»: учитывая желаемую дисперсионно-ковариационную матрицу которая по определению является положительно определенной, ее разложение Холецкого имеет вид: C = L L T ; L - нижняя треугольная матрица.$C$$C$$LL^T$$L$

Если вы теперь используете эту матрицу для проекции некоррелированного вектора случайных величин X , результирующая проекция Y = L X будет проекцией коррелированных случайных величин.$L$$X$$Y = LX$

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

+1 к @ user11852 и @ jem77bfp, это хорошие ответы. Позвольте мне подойти к этому с другой точки зрения, не потому, что я думаю, что это обязательно лучше на практике , а потому, что я думаю, что это поучительно. Вот несколько важных фактов, которые мы уже знаем:

1. представляет собой наклон линии регрессиикогда оба Х и Y являютсястандартизированы, то есть N ( 0 , 1 ) , $r$$X$$Y$$\mathcal N(0,1)$

2. - доля дисперсии в Y, относящаяся к дисперсии в X , $r^2$$Y$$X$

   
   
   (также из правил для отклонений ):

3. дисперсия случайной величины, умноженная на константу, представляет собой константу в квадрате, умноженную на исходную дисперсию:
   $$\text{Var}[aX]=a^2\text{Var}[X]$$
4. дисперсии сложения , т. е. дисперсия суммы двух случайных величин (при условии, что они независимы) является суммой двух дисперсий:
   $$\text{Var}[X+\varepsilon]=\text{Var}[X]+\text{Var}[\varepsilon]$$

Теперь мы можем объединить эти четыре факта, чтобы создать две стандартные нормальные переменные, популяции которых будут иметь заданную корреляцию (точнее, ρ ), хотя сгенерированные вами выборки будут иметь выборочные корреляции, которые различаются. Идея состоит в том, чтобы создать псевдослучайную переменную X , которая является стандартной нормалью, N ( 0 , 1 ) , а затем найти коэффициент a и дисперсию ошибки v e , такую, что Y ∼ N ( 0 , a 2 + v е ) , где$r$$\rho$$X$$\mathcal N(0,1)$$a$$v_e$$Y \sim\mathcal N(0,a^2+v_e)$ . (Обратите внимание, что | a | должно быть ≤ 1, чтобы это работало, и, кроме того, a = r .) Таким образом, вы начинаете с r , который хотите; это твой коэффициент, а . Затем вы вычисляете дисперсию ошибки, которая вам понадобится, это 1 - r 2 . (Если ваше программное обеспечение требует, чтобы вы использовали стандартное отклонение, возьмите квадратный корень из этого значения.) Наконец, для каждогосгенерированного вамипсевдослучайного значения x i сгенерируйте псевдослучайное значение ошибки, e $a^2+v_e=1$$|a|$ $\le 1$$a=r$$r$$a$$1-r^2$$x_i$$e_i$, с соответствующей дисперсией ошибки , и вычислить коррелированную псевдослучайную переменную, y i , путем умножения и сложения. $v_e$$y_i$

Если вы хотите сделать это в R, следующий код может работать для вас:

correlatedValue = function(x, r){
  r2 = r**2
  ve = 1-r2
  SD = sqrt(ve)
  e  = rnorm(length(x), mean=0, sd=SD)
  y  = r*x + e
  return(y)
}

set.seed(5)
x = rnorm(10000)
y = correlatedValue(x=x, r=.5)

cor(x,y)
[1] 0.4945964

(Изменить: я забыл упомянуть :) Как я уже описал, эта процедура дает вам две стандартные нормальные коррелированные переменные. Если вам не нужны стандартные нормали, но вы хотите, чтобы переменные имели определенные средние значения (не 0) и SD (не 1), вы можете преобразовать их, не влияя на корреляцию. Таким образом, вы должны вычесть наблюдаемое среднее значение, чтобы убедиться, что среднее значение равно , умножить переменную на нужный вам SD и затем добавить среднее значение, которое вы хотите. Если вы хотите, чтобы наблюдаемое среднее значение обычно колебалось вокруг желаемого среднего, вы бы вернули начальную разницу обратно. По сути, это преобразование z-счета в обратном направлении. Поскольку это линейное преобразование, преобразованная переменная будет иметь ту же корреляцию с другой переменной, что и раньше. $0$

Опять же, это, в простейшей форме, позволяет только генерировать пару коррелированных переменных (это можно увеличить, но очень быстро), и, конечно, это не самый удобный способ выполнить работу. В R вы хотели бы использовать ? Mvrnorm в пакете MASS , потому что это проще и потому что вы можете генерировать много переменных с заданной матрицей корреляции населения. Тем не менее, я думаю, что стоит пройти этот процесс, чтобы увидеть, как некоторые базовые принципы реализуются простым способом.

В общем, это не простая вещь, но я считаю, что есть пакеты для многофакторной генерации нормальной переменной (по крайней мере, в R, см. mvrnormВ MASSпакете), где вы просто вводите ковариационную матрицу и средний вектор.

Есть и еще один «конструктивный» подход. Допустим, мы хотим смоделировать случайный вектор и у нас есть его функция распределения F ( x 1 , x 2 ) . Первый шаг - получить функцию предельного распределения; т.е. интегрировать F по всем x 2 : F X 1 ( x 1 ) = ∫ ∞ - ∞ F ( x 1 , x 2 ) d x $(X_1,X_2)$$F(x_1,x_2)$$F$$x_2$ Затем мы находим F - 1 X 1 - обратную функцию от F X 1 - и включаем случайную величину ξ 1, которая равномерно распределена на интервале [ 0 , 1 ] . На этом шаге мы создаем первую координату х 1 = F - 1 X 1 ( £ , ) .

Теперь, так как мы получили одну координату, необходимо подключить его к исходной функции распределения , а затем получить условную функцию распределения с условием х 1 = х 1 : F ( х 2 | X 1 = х 1 ) = Р ( х 1 , х 2$F(x_1,x_2)$$x_1=\hat{x}_1$ гдеFХ1является функцией плотности вероятности предельногоХ1распределение; то естьF ' X 1 (x1)=fX1(x1).

$\xi_2$$[0,1]$$\xi_1$$F(x_2 | X_1=\hat{x}_1)$$\hat{x}_2=(F(x_2 | X_1=\hat{x}_1))^{-1}(\xi)$$\hat x_2$$F(\hat x_2 | X_1=\hat{x}_1) = \xi$

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

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

$\{x_i\}_{i=1}^N$$\{y_i\}_{i=1}^N$$C$

$c_{old}=corr(\{x_i\},\{y_i\})$

$n_1$$n_2: 1 \leq n_{1,2} \leq N$

$x_{n_1}$$x_{n_2}$

$c_{new}=corr( \{x_i\},\{y_i\})$

$|C-c_{new}| < |C-c_{old}|$

$|C-c| < \epsilon$

${x_i}$

Удачи!