Как генерировать коррелированые двоичные переменные

Вы описываете случайный процесс , и мне он кажется сложным ... если вы исключили требование двоичного (0,1) и вместо этого указали ожидаемое значение и дисперсию, это можно было бы описать как генератор белого шума, питающийся через однополюсный фильтр нижних частот, который, я думаю, даст вам характеристику α ^{| ij |} .

Это действительно может соответствовать планке для mathoverflow.net, в зависимости от того, как это сформулировано. Позвольте мне попытаться спросить ....

update: Я спрашивал на mathoverflow.net о случае α ^{| i-j |} . Но, возможно, там есть какие-то идеи, которые можно адаптировать к вашему случаю.

Выразите распределение x _i как линейную комбинацию некоторых независимых базисных распределений f _j : x _i = a _i1 f ₁ + a _i2 f ₂ +. .. . Ограничим f _j независимыми переменными, равномерно распределенными в 0..1 или в {0,1} (дискретных). Давайте теперь выразим все, что мы знаем, в матричной форме:

Let X be the vector (x1, x2, .., xn)
Let A be the matrix (a_ij) of dimension (k,n) (n rows, k columns)
Let F be the vector (f1, f2, .., fk) 
Let P be the vector (p1, p2, .., pn)
Let R be the matrix (E[x_i,x_j]) for i,j=1..n
Definition of the X distribution: X = A * F
Constraint on the mean of individual X variables: P = A * (1 ..k times.. 1)
Correlation constraint: AT*A = 3R or 2R in the discrete case (because E[x_i x_j] = 
  E[(a_i1*f_1 + a_i2*f_2 + ...)*(a_j1*f_1 + a_j2*f_2 + ...)] =
  E[sum over p,q: a_ip*f_p*a_jq*f_q] = (since for p/=q holds E[f_p*f_q]=0)
  E[sum over p: a_ip*a_jp*f_p^2] =
  sum over p: a_ip*a_jp*E[f_p^2] = (since E[f_p^2] = 1/3 or 1/2 for the discrete case)
  sum over p: 1/3 or 1/2*a_ip*a_jp
And the vector consisting of those sums over p: a_ip*a_jp is precisely AT*A.

Теперь вам нужно решить два уравнения:

AT*A      = 3R (or 2R in the discrete case)
A*(1...1) = P

Решение первого уравнения соответствует нахождению квадратного корня из матрицы 3R или 2R. См., Например, http://en.wikipedia.org/wiki/Cholesky_factorization и вообще http://en.wikipedia.org/wiki/Square_root_of_a_matrix . Что-то вроде. также следует сделать и со вторым :)

Я прошу окружающих математиков поправить меня, потому что я вполне мог смешать AT A с A AT или сделать что-то еще более неправильное.

Чтобы сгенерировать значение x _i в виде линейной смеси базисных распределений, используйте двухэтапный процесс: 1) используйте однородную случайную величину для выбора одного из базисных распределений, взвешенные с соответствующей вероятностью, 2) генерировать результат с использованием выбранного базисного распределения.

Быстрый поиск в RSeek показывает, что R имеет пакеты

• bindata
• binarySimCLF

для этого.

Решением грубой силы является выражение ограничений задачи в виде линейной программы с 2 ^ N переменными pr (w) , где w выходит за пределы все двоичные строки длины N . Во-первых, ограничение, что pr будет распределением вероятностей:

for all w: 0 <= pr(w) <= 1
sum_w pr(w) = 1

Во-вторых, ограничение, что математическое ожидание каждой переменной будет p :

for all i: sum_{w such that w[i] = 1} pr(w) = p

В-третьих, ковариационные ограничения:

for all i < j: sum_{w such that w[i] = w[j] = 1} pr(w) = const * |j - i|^alpha - p^2

Это очень медленно, но беглый поиск литературы не дал ничего лучшего. Если вы решите реализовать это, вот некоторые решатели LP с привязками Python: http://wiki.python.org/moin/NumericAndScientific/Libraries

Вот интуитивный / экспериментальный подход, который, похоже, работает.

Если b - двоичный ряд, m - среднее значение двоичной переменной, c - искомая корреляция, rand() генерирует U(0,1) r.v., и d - нужная вам корреляционная двоичная переменная:

d = if(rand() < c, b, if(rand() < m , 0, 1))

То есть если равномерная переменная меньше желаемой корреляции, d = b. В противном случае d = другое случайное двоичное число.

Я прогнал это 1000 раз для столбца из 2000 двоичных r.v. с m=.5 и c = .4 и c = .5. Среднее значение корреляции было точно таким, как указано, распределение оказалось нормальным. Для корреляции 0,4 среднеквадратичное отклонение корреляции составило 0,02.

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