# [R] Writing an R fn to produce random sequences of binaries XXXX

Duncan Murdoch murdoch.duncan at gmail.com
Tue Oct 15 14:17:14 CEST 2013

```On 13-10-15 7:55 AM, Dan Abner wrote:
> Hi all,
>
> I am attempting to write an R fn to produce a random sequence of 8 binaries
> (4 ordered pairs). I would like to parmeterize the fn so that it takes 2
> arguments: 1 that is an overall probability of a 1 vs. 0 and the other
> which controls the likelihood of 1s on the 1st vs. the 2nd element of the 4
> ordered pairs.
>
> Here are some examples:
>
>
>     *Voiced Probability Parm* *1st Element Probability Parm* *Average N of 1s
> * *Average N 1st Element 1s* * * *1* *1.5* *2* *2.5* *3* *3.5* *4*
> *4.5*
> 0.5 1 4 4 1 0 1 0 1 0 1 0  0.5 0 4 0 0 1 0 1 0 1 0 1  0.5 0.75 4 3 1 0 0 1 1
> 0 1 0  0.5 0.25 4 1 0 1 0 1 1 0 0 1  1 1 8 4 1 1 1 1 1 1 1 1  0 0 0 0 0 0 0
> 0 0 0 0 0  0.5 0.5 4 2 1 0 0 1 1 0 0 1
>
>
> Can anyone suggest an algorithm for this? I have a 1st draft of the fn
> definition, but my algorithm is not correct as I am not getting the average
> number of 1s (averaged over 10,000 sequences) correct.

A simple algorithm is to sample the 1st bit using rbinom(), then sample
the second bit, again using rbinom, but with probability conditional on
the value of the first bit.  But that requires a joint distribution, and
you haven't given us enough information to construct one.

For example, here's a model where bit 2 is either independent of bit 1
or forced to be equal to it, with theta being the probability of being
forced to be equal.  This probably isn't the model you want, but you
would just modify the bit2 line to give the correct conditional
distribution based on what you do want.

p <- 0.7     # marginal prob of a 1
theta <- 0.5 # ratio of times bit 2 is equal to bit 1, versus
# independent
bit1 <- rbinom(4, size=1, p)
bit2 <- rbinom(4, size=1, theta*bit1 + (1-theta)*p)
cbind(bit1, bit2)

Duncan Murdoch

>
> Thanks,
>
> Dan
>
> 	[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help