Re: generating non-IID random sequences
- To: mathgroup at smc.vnet.net
- Subject: [mg79203] Re: generating non-IID random sequences
- From: Mark Fisher <particlefilter at gmail.com>
- Date: Sat, 21 Jul 2007 04:21:02 -0400 (EDT)
- References: <f7kegd$58u$1@smc.vnet.net><f7poja$pac$1@smc.vnet.net>
On Jul 20, 3:35 am, Mark Fisher <particlefil... at gmail.com> wrote:
> On Jul 18, 3:12 am, Yaroslav Bulatov <yarosla... at gmail.com> 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 are 3 functions. The first is your version; the second is your
> version compiled; the third is a slightly different version also
> compiled.
>
> fun = Function[{n, p},
> NestList[
> RandomChoice[{p^# (1 - p)^(1 - #), p^(1 - #) (1 - p)^#} -> {1,
> 0}] &, RandomChoice[{0, 1}], n]];
>
> fun1 = Compile[{{n, _Integer}, {p, _Real}},
> NestList[
> RandomChoice[{p^# (1 - p)^(1 - #), p^(1 - #) (1 - p)^#} -> {1,
> 0}] &, RandomChoice[{0, 1}], n]];
>
> fun2 = Compile[{{n, _Integer}, {p, _Real}},
> NestList[
> If[# == 0, If[RandomReal[] < p, 0, 1],
> If[RandomReal[] < 1 - p, 0, 1]] &, RandomChoice[{0, 1}], n]];
>
> This gives the relative timings:
>
> First[#]/#& @ Table[Timing[i[10^5, .9]][[1]], {i, {fun, fun1, fun2}}]
>
> On my machine I get a factor of 5 speed up for fun1 and a factor of 25
> speed up for fun2 relative to fun.
>
> --Mark.
Well, I like both Peter's and Darren's solutions. Peter's can be
speeded up just a touch by compiling, but it beats mine even without
compiling. Darren's can be speeded up by using Peter's idea and
rearranging a bit. It ends up being about the same speed as Peter's.
markovBits1 =
Compile[{{n, _Integer}, {p, _Real}},
Module[{r =
Join[{RandomInteger[]}, RandomChoice[{p, 1 - p} -> {1, 0}, n]]},
Do[If[r[[i]] === 0, r[[i + 1]] = 1 - r[[i + 1]]], {i, n}];
r]]
--Mark