> -----Original Message-----
> From: Davidov, Ori [mailto:ori_davidov@merck.com]
> Sent: Thursday, 24 June 1999 4:10
> To: 's-news@wubios.wustl.edu'
> Subject: [S] Generating correlated binary random variables
>
>
>
> Hi,
>
> Does anyone out there have or knows of a S-Plus function that
> generates
> correlated binary rv's. Obviously there are many dependence
> structures and
> therefore many possible functions.
>
> Any information will be helpful to me at this point, thanks,
>
> Ori
Here's one idea.
Take X ~ Bernoulli(p), Y ~ Bernoulli(q), W ~ Bernoulli(r) and define
Z = (X|Y)&Z
Then Z ~ Bernoulli((p+(1-p)q)r) and
Cov(X,Z) = p(1-p)(1-q)r so that
Cor(X,Z) = sqrt(p(1-p)r) (1-q) / srqt[(p+(1-p)q) (1-pr-(1-p)qr)]
You can tune q and r to give desired values of E(Z) and Cor(X,Z). For
negative correlations use
Z = ((!X)|Y)&Z. The reason I used two auxiliary random variables Y and W was
to get a truth table for X and Z with all cells having non-zero probability
(see tables in example below).
Examples follow:
> p<-q<-r<-.5
> x<-rbinom(1000,1,p)
> y<-rbinom(1000,1,q)
> w<-rbinom(1000,1,r)
> z<-(x|y)&w
> apply(cbind(x,y,w,z),2,mean)
x y w z
0.511 0.473 0.48 0.356
> (p+(1-p)*q)*r # E(z)
[1] 0.375
> round(var(cbind(x,y,w,z)),2)
x y w z
x 0.25 -0.01 0.00 0.06
y -0.01 0.25 -0.01 0.05
w 0.00 -0.01 0.25 0.19
z 0.06 0.05 0.19 0.23
> p*(1-p)*(1-q)*r # Cov(x,z)
[1] 0.0625
> round(cor(cbind(x,y,w,z)),2)
x y w z
x 1.00 -0.03 -0.01 0.26
y -0.03 1.00 -0.04 0.20
w -0.01 -0.04 1.00 0.77
z 0.26 0.20 0.77 1.00
> sqrt(p*(1-p)*r)* (1-q) / sqrt((p+(1-p)*q)* (1-p*r-(1-p)*q*r)) # Cor(x,z)
[1] 0.2581989
> table(x,x|y)
FALSE TRUE
0 250 239
1 0 511
> table(x,x&y)
FALSE TRUE
0 489 0
1 277 234
> table(x,z)
FALSE TRUE
0 376 113
1 268 243
Nick Ellis
CSIRO Marine Research mailto:Nick.Ellis@marine.csiro.au
PO Box 120 ph +61 (07) 3826 7260
Cleveland QLD 4163 fax +61 (07) 3826 7222
Australia http://www.marine.csiro.au
-----------------------------------------------------------------------
This message was distributed by s-news@wubios.wustl.edu. To unsubscribe
send e-mail to s-news-request@wubios.wustl.edu with the BODY of the
message: unsubscribe s-news
|