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>
- generating non-IID random sequences