Let list be a list of lists, the elements of which are random integers in the interval $[1,c]$. Each element $i$ has a probability prob to be changed to another integer, selected randomly in $[1,c] \setminus \{v_i\}$ where $v_i$ is the value of the element $i$. No particular size is requested for either list or its sublists, but each sublist has the same size.

As an example, let's suppose that (a) list is made of 3 sublists, (b) each sublist contains 2 elements, (c) c equals to 3, and (d) prob is 0.25.

A random generation would give for instance list = {{3, 1}, {2, 1}, {3, 3}}. After "running" the probability through the elements of the sublists, we would have to change for instance the elements at positions {1, 2} and {3, 1}. So, 1 will be randomly changed to either 2 or 3, and 3 to either 1 or 2. A possible outcome of list = {{3, 1}, {2, 1}, {3, 3}} would then be {{3, 2}, {2, 1}, {2, 3}}.

I am interested in an efficient way of coding this problem, with the specifications given in the first paragraph.

I tried a method with three different versions, shown below. The respective computing times are given for a problem where (a) list is made of 1000 sublists, (b) each sublist contains 30 elements, (c) c equals to 3, and (d) prob is 0.001.

list = RandomInteger[{1, 3}, {1000, 30}]; 

Version 1

ChangeV1[prob_, c_, sizelist_, sizesublist_][list_] := 
 Block[{elemstochange = RandomVariate[BernoulliDistribution[prob], {sizelist, sizesublist}]},
 ReplacePart[list, (# -> RandomChoice[Drop[Range[c], {Extract[list, #]}]]) 
                     & /@ Position[elemstochange, 1]]];


ChangeV1[0.001, 3, 1000, 30][list]; // AbsoluteTiming
(* {0.00297177, Null} *)

Version 2

ChangeV2[prob_, c_, sizelist_, sizesublist_][list_] := 
 Block[{elemstochange = RandomVariate[BernoulliDistribution[prob], {sizelist, sizesublist}]},
 MapAt[RandomChoice[Drop[Range[c], {#}]] &, list, Position[elemstochange, 1]]];


ChangeV2[0.001, 3, 1000, 30][list]; // AbsoluteTiming
(* {0.00313513, Null} *)

Version 3

ChangeV3[prob_, c_, sizelist_, sizesublist_][list_] := 
 Block[{elemstochange = RandomVariate[BernoulliDistribution[prob], {sizelist, sizesublist}],
        temp = list},
 (temp[[Sequence @@ #]] = 
  RandomChoice[Drop[Range[c], {Extract[list, #]}]]) & /@ Position[elemstochange, 1];


ChangeV3[0.001, 3, 1000, 30][list]; // AbsoluteTiming
(* {0.0028965, Null} *)

I would be very interested in reading other methods that would improve the computing time for Change.

  • $\begingroup$ Is the end result all you're after (hence use of e.g. Nest )? If so, seems like a lot of machinations that could be done simply. How big might c be? "Small" (e.g. <thousands) or does it have a large upper bound? What of the probabilities? Will they be "small" as in your test cases, or do they cover the [0,1] space completely? It would be useful to know more about "...main routine I am working on, there are other functions involved that make modifications..." - what is th purpose of this and the overall goal for the "main" routine - there might be a much more efficient way to get there. $\endgroup$
    – ciao
    Aug 25, 2015 at 5:21
  • $\begingroup$ I reformulated the original question to make things clearer. $\endgroup$
    – user31159
    Aug 25, 2015 at 21:41

2 Answers 2


Here's a single-shot adaptation of my earlier (deleted since question changed) ideas. A bit over twice as fast as those for single-shot cases only, e.g. on c=5, 5000 sublists, 30 sublist length, p=0.001 done once.

mutateList5[list_, c_, p_] := 
 Module[{fl = Flatten@list, cr = Range@c, cl, lfl = Times @@ Dimensions@list, chg},
  cl = Delete[cr, #] & /@ cr;
  chg = RandomSample[Range@lfl, 
    RandomVariate[BinomialDistribution[lfl, p]]];
  fl[[chg]] = Extract[cl, Transpose[{fl[[chg]], RandomInteger[{1, c - 1}, Length@chg]}]];
  ArrayReshape[fl, Dimensions@list]];

This can be quite a bit faster than above (as in using 2000 sublists of size 30 with c of 10 and p of 0.1 it was over 50X faster than ChangeV3 and over 40X faster than ChangeV1) :

mutateList7[list_, c_, p_, method_: 1] := 
  Block[{th = 50, dl = Dimensions@list, lfl, fl, cr, chg, cn},
   lfl = Times @@ dl;
   cn = Switch[method, 1, RandomVariate[BinomialDistribution[lfl, p]],
                       2, Tr@RandomVariate[BernoulliDistribution[p], lfl],
                       3, RandomVariate[PoissonDistribution[lfl*p]]];
   Which[cn == 0, list,
    cn <= th,
    chg = RandomSample[Range@lfl, cn];
    chg = Transpose[{Ceiling[chg/dl[[2]]], Mod[chg, dl[[2]], 1]}];
    cr = Range@(c - 1);
    ReplacePart[list, Thread[chg -> Mod[Extract[list, chg] + RandomChoice[cr, cn], c, 1]]],
    fl = Flatten@list;
    cr = Range@(c - 1);
    chg = RandomSample[Range@lfl, cn];
    fl[[chg]] = Mod[fl[[chg]] + RandomChoice[cr, cn], c, 1];
    ArrayReshape[fl, Dimensions@list]]];

Uses heuristics to optimize strategy for changes. This also uses an optional last argument that uses differing methods for change determination. I've found some... interesting behavior of Binomial RV generation, where it was much slower for your test case than made sense, and summing Bernoulli was much faster, yet on larger cases it was fast as expected. I'm investigating this, I'd venture it's some auto-compilation or internal algorithm switch depending on parameters.

In any case, try them out. Method 1 (default) uses Binomial, 2 uses sum of Bernoulli, and method 3 uses Poisson RV, which can be very fast and is appropriate when number of elements * probability <= 10 to 20 or so.

As an aside, you've still not clarified what happens to parameters between calls for a given "run" of your "thing", i.e., does c change? Does p change? Is it called with differing contents and perhaps dimensions, but p or c, or both, are constant? If any of these are true, significant speed gains can be had...

Update: Per comments, here's a potentially much more efficient means to the end:

mutateBldr[c_, p_, sl_, sls_, mr_] := 
  Module[{rv, rvf, rvc = 1, ts, cr, fn, th = 100},
   rv = RandomVariate[BinomialDistribution[sl sls, p], mr];
   ts = Tuples[{Range@sl, Range@sls}];
   rvf = Range@(sl*sls);
   cr = Range@(c - 1);
   fn = If[rv[[rvc]] == 0, rvc++; #, Block[{t, m},
       If[rv[[rvc]] > th, m = Flatten@#;
        t = RandomSample[rvf, rv[[rvc]]];
        m[[t]] = Mod[m[[t]] + RandomChoice[cr, rv[[rvc++]]], c, 1];
        Partition[m, sls],
        t = RandomSample[ts, rv[[rvc]]];
        ReplacePart[#, Thread[t -> Mod[Extract[#, t] + RandomChoice[cr, rv[[rvc++]]], c, 1]]]]]] &;

Use: Prime it with the specific c, # sublists, sublist size, and probability, along with expected runs (so if you're going to call it say 1000 times, use 1500 or so to give wiggle-room):

mymute = mutateBldr[c, prob, subLists, subListSz, 10000];

Then, just call it on your evolving lists:


I think you'll find this significantly outperforms prior methods...

  • $\begingroup$ That's great, thanks a lot @ciao. The method with Mod is nice and the use of a Poisson distribution is a good idea. About your questions, for a given run of the "main" function, c, p and the dimensions of list are constant. Only the contents of the sublists within list can change. (One way they change is from mutations (with mutateList), the other is from crossovers – I am considering a genetic algorithm). $\endgroup$
    – user31159
    Aug 26, 2015 at 14:12
  • $\begingroup$ @Xavier Cool. So, if you need that last sliver of speed, you could build a function that takes p c and dims, precomputes the quantiles and other giblets and returns a function that you use for runs for those parameters. I'd venture you could double performance or so... $\endgroup$
    – ciao
    Aug 26, 2015 at 19:28
  • $\begingroup$ @Xavier - just cobbled up a quick-n-dirty of a pre-computed version, turns out it's over an order of magnitude faster... worth pursuing for you prehaps. $\endgroup$
    – ciao
    Aug 26, 2015 at 23:11
  • $\begingroup$ Ah that's very nice @ciao. Thanks for the tip, I'll add that. $\endgroup$
    – user31159
    Aug 26, 2015 at 23:12
  • $\begingroup$ @Xavier: Added the pre-computed version, give it a whirl and report back... $\endgroup$
    – ciao
    Aug 27, 2015 at 4:36

In the narrowly defined interpretation of the question, one can get an incremental speed gain by flattening the lists, and reshaping the result to the desired dimensions afterwards:

ChangeV4[prob_, c_, sizelist_, sizesublist_][list_] := 
  Block[{elemstochange = 
      sizelist sizesublist],
    temp = Flatten@list,
    choices = Table[Drop[Range[c], {i}], {i, c}]
   (temp[[#]] = RandomChoice[choices[[temp[[#]]]]]) & /@ 
    Flatten@Position[elemstochange, 1];
   Partition[temp, sizesublist]];

Here, I also tried to get another small gain in efficiency by taking the operation Drop[Range[c],...] out of the "loop" over the elements to change.

In my test with list as given in the question, this performed faster than all three original versions.

  • $\begingroup$ Thanks for your answer @Jens and the ideas with Flatten and Drop. $\endgroup$
    – user31159
    Aug 26, 2015 at 13:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy