模拟重复二进制测量的数据

Question

我可以生成一个二进制变量y如下：

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

如何生成 n 个变量 y_1, y_2, ..., y_n，相关性为 rho？

Answer 1

相关性是成对衡量，所以我假设当你谈论二进制（伯努利）值 Y₁,...,Y_n 具有 rho 的相关性，您将它们视为时间序列 Y_i: i = 1,...,n，伯努利值具有共同均值 p、方差 p*(1-p) 和 rho.

的滞后 1 相关性

我能够使用相关性和条件概率的定义来计算它。鉴于这是一堆乏味的代数，而且 Whosebug 没有优雅地进行数学计算，我直接跳到结果，用伪代码表示：

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 它只生成 Bernoulli(p)，而不管先前的值。正如您在问题中已经指出的那样，伯努利 RV 是 n = 1.

的二项式

这适用于所有 0 <= rho, p <= 1。对于负相关，p和rho的相对大小有限制，因此伯努利的参数总是在0和1之间。

您可以分析检查条件概率以确认正确性。我不使用 Stata，但我在 JMP statistical software 中对其进行了非常彻底的测试，它非常有效。

实施（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 [分别在左下方和右下方]。这些估计值与具有 rho = 0.5 的 Bernoulli(0.7) 完全匹配，如演示程序中指定的那样。

如果目标是生成具有指定相关性的 (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)

Answer 2

这是@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。

Answer 3

这是@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          |                  |

模拟重复二进制测量的数据

Simulate data for repeated binary measures

simulation

variables

stata