Services & ResourcesWolfram ForumsMathGroup Archive

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

Re: generating non-IID random sequences

• To: mathgroup at smc.vnet.net
• Subject: [mg79230] Re: generating non-IID random sequences
• From: Yaroslav Bulatov <yaroslavvb at gmail.com>
• Date: Sat, 21 Jul 2007 04:35:01 -0400 (EDT)
• References: <f7kegd$58u$1@smc.vnet.net><f7poja$pac$1@smc.vnet.net>

The Compiled one-liner (fun2) seems to be the fastest solution.

I also compared it with a simple java solution, and it was about 2.7
times faster (19 seconds vs 7 seconds). Using Java solution through
JLink seems to take about the same time as running stand-alone.
Running JVM with some optimization flags (server,compileThreshold)
took the time to 6 seconds. A simple C solution took 2 seconds, 1.5
seconds with gcc -O9 optimization.

The bottom line seems to be that efficient Mathematica solution can be
25 times faster than inefficient one, porting to C can give 10x speed-
up, and Java is somewhere in between.


fun2 = Compile[{{n, _Integer}, {p, _Real}},
    NestList[
     If[# == 0, If[RandomReal[] < p, 0, 1],
       If[RandomReal[] < 1 - p, 0, 1]] &, RandomChoice[{0, 1}], n]];
Timing[fun2[4 10^7, .9]][[1]]

public class hmm {
    public static int[] sample(double p, int length) {
	int[] result = new int[length];
	result[0]=(Math.random()>.5)?0:1;
	for (int i = 1; i < length; ++i)
	    result[i]=(Math.random()>p)?(1-result[i-1]):result[i-1];
	return result;
    }

    public static void main(String[] args) {
        long start = System.currentTimeMillis();
        sample(.9,4*10000000);
        System.out.println((System.currentTimeMillis()-start)/1000.);
    }
}


#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define n 40000000

main() {
  int* result = malloc(4*n);
  double start = clock();
  int i;
  result[0] = 0;
  for (i = 1; i<n;++i)
    result[i]=(rand()/((double)RAND_MAX + 1)>.9)?(1-result[i-1]):
(result[i-1]);
  printf("%d\n",result[i-1]);
  printf("%f\n",(clock()-start)/CLOCKS_PER_SEC);
  return 0;
}


On Jul 20, 12: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.

• Follow-Ups:
  • Re: generating non-IID random sequences
    • From: Mark Fisher <particlefilter@gmail.com>
  • Re: Re: generating non-IID random sequences
    • From: Carl Woll <carlw@wolfram.com>