Корреляция является попарной мерой, поэтому я предполагаю, что когда вы говорите о двоичных (Бернулли) значениях, Y 1 , ..., Y n , имеющих корреляцию rho
Вы просматриваете их как временной ряд Y i : i = 1, ..., n, значений Бернулли, имеющих общее среднее p
, дисперсию p*(1-p)
и лаг 1корреляция rho
.
Мне удалось решить это, используя определение корреляции и условной вероятности.Учитывая то, что это была куча утомительной алгебры, а stackoverflow не справляется с математикой изящно, я сразу перехожу к результату, выраженному в псевдокоде:
if Y[i] == 1:
generate Y[i+1] as Bernoulli(p + rho * (1 - p))
else:
generate Y[i+1] as Bernoulli(p * (1 - rho))
В качестве проверки работоспособности вы можете увидеть, что если rho = 0
он просто генерирует Бернулли (р) независимо от предыдущего значения.Как вы уже отметили в своем вопросе, RV Бернулли являются биномами с n = 1
.
. Это работает для всех 0 <= rho, p <= 1
.Для отрицательных корреляций существуют ограничения на относительные величины p
и rho
, так что параметры Бернулли всегда находятся в диапазоне от 0 до 1.
Вы можете аналитически проверить условные вероятности, чтобы подтвердить правильность.Я не использую Stata, но я довольно тщательно проверил это в статистическом программном обеспечении JMP , и оно работает как шарм.
IMPLEMENTATION (Python)
import random
def Bernoulli(p):
return 1 if random.random() <= p else 0 # yields 1 w/ prob p, 0 otherwise
N = 100000
p = 0.7
rho = 0.5
last_y = Bernoulli(p)
for _ in range(N):
if last_y == 1:
last_y = Bernoulli(p + rho * (1 - p))
else:
last_y = Bernoulli(p * (1 - rho))
print(last_y)
Я запустил это и перенаправил результаты в файл, затем импортировал файл в JMP.Анализируя его как полученный временной ряд:
Среднее значение выборки составляло 0,69834 со стандартным отклонением 0,4589785 [верхняя правая часть рисунка].Оценки lag-1 для автокорреляции и частичной корреляции составляют 0,5011 [слева внизу и справа соответственно].Все эти оценочные значения являются отличными совпадениями с Бернулли (0,7) с rho = 0.5
, как указано в демонстрационной программе.
Если вместо этого целью является получение (X, Y) пар с указаннымкорреляция, измените цикл на:
for _ in range(N):
x = Bernoulli(p)
if x == 1:
y = Bernoulli(p + rho * (1 - p))
else:
y = Bernoulli(p * (1 - rho))
print(x, y)