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

Мне нужно создать серию из N случайных двоичных переменных с заданной корреляционной функцией. Пусть x = {x _i} - это серия двоичных переменных (с использованием значения 0 или 1, я от 1 до N). Маргинальной вероятности задается Pr (x _i= 1) = p, а переменные должны быть скоррелированы следующим образом:

Corr [x _i x _j] = const & times; | я & minus; j | ^{& minus; & alpha;} (для я!= j)

где & alpha; - положительное число.

Если это проще, рассмотрим корреляционную функцию:

Corr [x _i x _j] = (| я & minus; j | +1) ^{& minus; & alpha;}

Существенная часть заключается в том, что я хочу исследовать поведение, когда корреляционная функция идет как степенной закон. (не & alpha; ^{| я & minus; j |})

Можно ли создать такую серию, предпочтительно в Python?

Ответ 1

Спасибо за все ваши материалы. Я нашел ответ на свой вопрос в симпатичной маленькой статье Chul Gyu Park и др., Поэтому, если кто-то столкнется с той же проблемой, посмотрите:

"Простой способ генерации коррелированных двоичных переменных" (jstor.org.stable/2684925)

для простого алгоритма. Алгоритм работает, если все элементы в корреляционной матрице положительны и для общего маргинального распределения Pr (x_i) = p_i.

Ответ 2

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

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

update: я сделал запрос на mathoverflow.net для случая & alpha; ^{| i-j |}. Но, возможно, есть некоторые идеи, которые могут быть адаптированы к вашему делу.

Ответ 3

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

чтобы сделать это.

Ответ 4

Выразите распределение 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. Что-то должно быть сделано и во втором:)

Я прошу математиков исправить меня, потому что я вполне мог бы смешать ATA с AAT или сделать что-то еще более неправильное.

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

Ответ 5

Решение грубой силы должно выражать ограничения задачи как линейной программы с переменными 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

Ответ 6

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

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

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

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

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

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