Имитация данных для повторных бинарных измерений - PullRequest
0 голосов
/ 27 ноября 2018

Я могу сгенерировать двоичную переменную y следующим образом:

clear
set more off
gen y =.
replace y = rbinomial(1, .5)

Как я могу сгенерировать n переменных y_1, y_2, ..., y_n с корреляцией rho?

Ответы [ 3 ]

0 голосов
/ 30 ноября 2018

Это решение @ pjs в Stata для генерации пар переменных:

clear
set obs 100
set seed 12345

generate x = rbinomial(1, 0.7)

generate y = rbinomial(1, 0.7 + 0.2 * (1 - 0.7)) if x == 1
replace y  = rbinomial(1, 0.7 * (1 - 0.2)) if x != 1

summarize x y

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
           x |        100         .72    .4512609          0          1
           y |        100         .67    .4725816          0          1

correlate x y
(obs=100)

             |        x        y
-------------+------------------
           x |   1.0000
           y |   0.1781   1.0000

И симуляции:

set seed 12345

tempname sim1
tempfile mcresults

postfile `sim1' mu_x mu_y rho using `mcresults', replace

forvalues i = 1 / 100000 {
    quietly {
        clear
        set obs 100

        generate x = rbinomial(1, 0.7)
        generate y = rbinomial(1, 0.7 + 0.2 * (1 - 0.7)) if x == 1
        replace  y = rbinomial(1, 0.7 * (1 - 0.2)) if x != 1

        summarize x, meanonly
        scalar mean_x = r(mean)

        summarize y, meanonly
        scalar mean_y = r(mean) 

        corr x y
        scalar rho = r(rho)

        post `sim1'  (mean_x) (mean_y) (rho)
    }
}

postclose `sim1'

use `mcresults', clear

summarize *

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
        mu_x |    100,000    .7000379    .0459078        .47        .89
        mu_y |    100,000    .6999094    .0456385        .49        .88
         rho |    100,000    .1993097    .1042207  -.2578483   .6294388

Обратите внимание, что в этом примере я использую p = 0.7 и rho = 0.2.

0 голосов
/ 09 декабря 2018

Это решение @ pjs в Stata для генерации временных рядов :

clear
set seed 12345
set obs 1

local p = 0.7
local rho = 0.5

generate y = runiform()
if y <= `p' replace y = 1
else replace y = 0

forvalues i = 1 / 99999 {
    set obs `= _N + 1'
    local rnd = runiform()
    if y[`i'] == 1 {
        if `rnd' <= `p' + `rho' * (1 - `p') replace y = 1 in `= `i' + 1'
        else replace y = 0 in `= `i' + 1' 
    }
    else {
        if `rnd' <= `p' * (1 - `rho') replace y = 1 in `= `i' + 1'
        else replace y = 0 in `= `i' + 1' 
    }
}

Результаты:

summarize y

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
           y |    100,000      .70078    .4579186          0          1


generate id = _n
tsset id
corrgram y, lags(5)

                                          -1       0       1 -1       0       1
 LAG       AC       PAC      Q     Prob>Q  [Autocorrelation]  [Partial Autocor]
-------------------------------------------------------------------------------
1        0.5036   0.5036    25366  0.0000          |----              |----    
2        0.2567   0.0041    31955  0.0000          |--                |        
3        0.1273  -0.0047    33576  0.0000          |-                 |        
4        0.0572  -0.0080    33903  0.0000          |                  |        
5        0.0277   0.0032    33980  0.0000          |                  |        
0 голосов
/ 29 ноября 2018

Корреляция является попарной мерой, поэтому я предполагаю, что когда вы говорите о двоичных (Бернулли) значениях, 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.Анализируя его как полученный временной ряд:

See description below

Среднее значение выборки составляло 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)
...