Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

Re: generating non-IID random sequences

• To: mathgroup at smc.vnet.net
• Subject: [mg79124] Re: [mg79094] generating non-IID random sequences
• From: Darren Glosemeyer <darreng at wolfram.com>
• Date: Thu, 19 Jul 2007 03:30:29 -0400 (EDT)
• References: <200707180701.DAA04619@smc.vnet.net>

Yaroslav Bulatov wrote:
> I'm looking for a fast way to sample from a Markov-1 sequence of
> random bits. The method below is 3600 times slower than built-in
> RandomInteger function, can it be made much faster?
>
> p = 0.9; (* the probability of encountering 00 or 11 *)
> f = RandomChoice[{p^# (1 - p)^(1 - #), p^(1 - #) (1 - p)^#} -> {1, 0}]
> &;
> NestList[f, RandomChoice[{0, 1}], 100000]
>
>   

Here's some code that is roughly 25-30 times faster than the code sent 
for a list of 10^6 bits. The code starts with a list containing a 0 or 1 
(generated with equal probability) followed by n 0 or 1 values generated 
with 1 having probability p. Those values are then looped through, 
iteratively replacing any value that follows a 0 with the other value 
(the probabilities flip-flop for values following a 0, so we can just 
switch the value). In effect this trades the more costly n calls to 
RandomChoice for a single call to RandomChoice and an update loop 
through the list, and uses compilation to speed up the pass through the 
loop.


In[1]:= markovBits[p_, n_] := Block[{not},
          Compile[{{rr, _Integer, 1}},
            not[1] = 0;
            not[0] = 1;
            Module[{rrnew = rr},
             Do[If[rrnew[[i]] === 0,
               rrnew[[i + 1]] = not[rrnew[[i + 1]]]], {i,
               Length[rrnew] - 1}];
             rrnew]
            ][Join[{RandomInteger[]}, RandomChoice[{p, 1 - p} -> {1, 0}, 
n]]]]

In[2]:= markovBits[.9, 10^6]; // Timing

Out[2]= {0.89, Null}


The following counts the values of absolute differences of adjacent 
values to give a check of the adjacency probabilities, so the 0s are the 
00 and 11 cases.

In[3]:= Tally[Abs[Most[#] - Rest[#]]] &[markovBits[.9, 10^6]]

Out[3]= {{0, 900210}, {1, 99790}}


Darren Glosemeyer
Wolfram Research

• References:
  • generating non-IID random sequences
    • From: Yaroslav Bulatov <yaroslavvb@gmail.com>