Efficiently Generating Lists in a Mathematica Way

Question

I'm struggling to move away from a procedural coding background and learn good Mathematica functionality. Most commonly, my issue is with generating lists where a little computation is required to determine the entry in each list. My understanding is that the Table function should be used in creating lists, but I struggle to make it work in these kind of situations. To give in example, suppose I am drawing random numbers from [0,10] and I wanted to estimate the expected value of how many draws are required to get a number > 5. One way would be to do a bunch of trials and for each trial see how many numbers I had to draw to get one > 5. Given the number of trials, numRums, I'd like to store this in a list. I can do this with a for loop:

HowManyTimesList = {};
numRuns = 10;
For[i = 1, i <= numRuns, i++,
 iter = 1;
 While[RandomReal[{0, 10}] < 5, iter++]
  AppendTo[HowManyTimesList, iter];
 ]
N[Mean[HowManyTimesList]]

But this seems ugly and un-Mathematica-esque. My struggle with using a Table is that I need to loop through doing draws until I have a number > 5, but I then need to return the number of times. I can't fit that into the Table function because it seems to only take single expression as the first argument.

I imagine there might be clever ways to organize this particular example so that it returns the desired expression in a single line, and then can be squeezed into a Table. But it seems like there should be a nice clean way to make lists where you don't always have to squeeze things into a single returnable line. So how do you approach such situations in a proper Mathematica way?

Answer 1

To answer your question about using Table (or Do) instead of For, here's a version of your code using Do:

SeedRandom[1];
HowManyTimesList = {};
numRuns = 10;
Do[
    iter=1; While[RandomReal[10]<5, iter++]; AppendTo[HowManyTimesList, iter],
    {numRuns}
]

HowManyTimesList

{1, 2, 4, 3, 2, 3, 1, 1, 1, 3}

If you use Table instead, you don't need to use AppendTo + HowManyTimesList:

SeedRandom[1];
numRuns = 10;
Table[
    iter=1; While[RandomReal[10]<5, iter++]; iter,
    {numRuns}
]

{1, 2, 4, 3, 2, 3, 1, 1, 1, 3}

Answer 2

Typically, to construct a list of results until some condition is (not) achieved, I would use NestWhileList.

This takes a function and repeatedly applies it to the previous output until the output does not satisfy a given condition. Something like:

NestWhileList[RandomReal[{0, 10}] &, RandomReal[{0, 10}], # <= 5 &]

Would output, perhaps:

{4.79448, 2.24396, 5.95265}

To explain each part of this to try and make it as clear as possible:

RandomReal[{0, 10}]& is an anonymous function which receives and ignores one input and returns the result of RandomReal[{0, 10}], the second RandomReal[{0, 10}] is the input value to NestWhileList, and # <= 5 & is a function which returns True if the condition you are looking for (# > 5 &) has not yet been met.

If you'd like to avoid just shoving this whole thing into a table, you could save it as a delayed expression using :=:

expr := NestWhileList[RandomReal[{0, 10}] &, RandomReal[{0, 10}], # <= 5 &]

Whenever expr is evaluated, then, it will re-run this process. Now to evaluate this some number of times, it's as simple as:

numRuns = 10;
Table[Length[expr], {numRuns}]

I approach it this way because it naturally lends itself to any other statistical values you might want given these exact lists as well. Say you want to find the average of the average values of the list members, you could do something like:

Mean[Table[Mean[expr], {numRuns}]]

Without any excessive repetition of code.

For other applications, I would also recommend looking at NestList, Fold, Inner, Outer, and Sow and Reap, each of which have lots of practical and very general uses in list manipulation.

Answer 3

numRuns = 10;

Two alternative approaches:

1. Use RandomVariate to generate the numRuns random numbers from GeometricDistribution[.5]:

1 + RandomVariate[GeometricDistribution[.5], numRuns]

{4, 2, 4, 1, 1, 1, 1, 3, 2, 4}

2. Generate numRuns random numbers at a time until we collect a list of runs with length numRuns:

Module[{lst = UnitStep[RandomReal[10, numRuns] - 5]}, 
 While[Total[lst = Join[lst, UnitStep[RandomReal[10, numRuns] - 5]]] < numRuns]; 
 Take[Differences[Join[{0}, Random`Private`PositionsOf[lst, 1]]], numRuns]]

{2, 1, 2, 3, 3, 2, 2, 3, 1, 3}

Timings:

For large numRuns both approaches are faster than the approaches generating a single random number at a time:

numRuns = 10000;

SeedRandom[1]
(res1 = 1 + RandomVariate[GeometricDistribution[.5], numRuns];) // AbsoluteTiming // First

0.0003668

SeedRandom[1];

(res2 = Module[{lst = UnitStep[RandomReal[10, numRuns] - 5]}, 
      While[Total[lst = Join[lst, UnitStep[RandomReal[10, numRuns] - 5]]] < numRuns]; 
      Take[Differences[Join[{0}, Random`Private`PositionsOf[lst, 1]]], numRuns]];) //
   AbsoluteTiming // First

0.0011162

compare with the methods from Carl's and eyorble's answers:

SeedRandom[1];
(res3 = Table[iter = 1; While[RandomReal[10] < 5, iter++]; iter, {numRuns}];) // 
  AbsoluteTiming // First

0.0230272

SeedRandom[1];
(res4 = Table[Length[expr], {numRuns}];) // AbsoluteTiming // First

0.041832

Some pictures:

ClearAll[pdf1, pdf2, pdf3, pdf4]
pdf1[x_?NumericQ] := PDF[EmpiricalDistribution@res1, x]
pdf2[x_?NumericQ] := PDF[EmpiricalDistribution@res2, x]
pdf3[x_?NumericQ] := PDF[EmpiricalDistribution@res3, x]
pdf4[x_?NumericQ] := PDF[EmpiricalDistribution@res4, x]

Row[{BarChart[Through[{pdf1, pdf2, pdf3, pdf4}@#] & /@ Range[0, 10], 
    ChartLayout -> "Grouped", Frame -> True,  ChartLabels -> {Range[0, 10], None}, 
    ChartStyle -> 97, ImageSize -> Medium],
  BoxWhiskerChart[{res1, res2, res3, res4}, ChartStyle -> 97, 
    ChartLegends -> {"res1", "res2", "res3", "res4"}, 
    ChartLabels -> {"res1", "res2", "res3", "res4"}, ImageSize -> Medium]}, Spacer[5]]