# Random Sample from Very Large Table

I'm currently trying to create some test data to check a numerical implementation against Mathematica's implementation of some functions.

In order to do this, I've been generating a broad range of arguments with a (Parallel)Table which works fine for a function with a few arguments, but I have one function I want to test which takes 4 integer arguments and 10 floating point arguments:

f[n1, n2, n3, n4, x0, x1, x2, x3, x4, x5, x6, x7, x8, x9]


The thing is creating a broad range of values use a table results in a simply enormous amount of data which cannot be held in memory. For example even:

realSample[n_] := Join[-10^Subdivide[-10, 10, n], 10^Subdivide[-10, 10, n], Subdivide[-10, 10, n]];
Table[{n1, n2, n3, n4, a, b, c, d, e, f, w, x, y, z},
{n1, 0, 1}, {n2, 0, 1}, {n3, 0, 1}, {n4, 0, 1},
{x0, realSample[1]}, ...
]


results in something that is just unmanageable even for n = 1 simply because of the sheer number of arguments.

So instead, I'd like to take a random sample of the whole table (say a million); however, RandomSample doesn't stream the table and the same problem as before still comes up.

Short of implementing a reservoir sampling, is there anything in Mathematica that allows me to take a random sample of the full table without loading the whole table into memory?

• May be apply RandomSample[] on each of the arguments separately and then build the full argument? RandomSample[#,1]&/@{{0,1}, {0,1}, ...}. Jun 6 '20 at 5:42
• Although this can work, it's not quite equivalent to random sampling the whole table Jun 6 '20 at 5:50
• RandomSample[] is surely the correct thing to use, but perhaps use it as Join[RandomSample[{0, 1}, 4], RandomSample[realSample[n], 10]]? Jun 6 '20 at 6:27
• Why don't you just RandomSample the positions of entries in the table and then create/load the entries of the table only for the sampled positions? Jun 6 '20 at 6:27
• @HenrikSchumacher That actually sounds very reasonable and I hadn't thought of that! Jun 6 '20 at 7:40

So one solution alluded to in the question is to use reservoir sampling which only holds in memory the desired final number of items while iterating through the table:

Attributes[ReservoirSample] = {HoldRest};
ReservoirSample[n_, arg_, iter__] := Block[{
sample = {}, p = 0
},

Do[
p += 1;
If[Length[sample] < n,

(* Fill the reservoir to begin with *)
AppendTo[sample, arg]
,

(* Otherwise we add the next argument with decreasing probability *)
If[RandomInteger[{1, p}] <= n,
sample[[RandomInteger[{1, n}]]] = arg;
];
];
,
iter
];

sample
];


Checking the output in a very simple case does seem to produce a flat distribution as desired:

Join @@ Table[ReservoirSample[
10,
x,
{x, 10},
{y, 10},
{z, 10}
], 100] // Histogram


I have tried implementing a parallel reservoir sampling, but due to the way Mathematica subdivides the iteration between the kernels, I'm consistently getting samples which are very much non-uniform.

Following the suggestion of @HenrikSchumacher, here is another solution which creates the entries based on randomly selecting a value for each iteration variable at each step.

This is much faster than reservoir sampling in my other answer as it does not need to iterate through every possibility. The main drawback however is that this can produce duplicates, especially if the number of samples requested isn't much smaller than the total number of possibilities. The reservoir sampling is guaranteed to not do that.

Attributes[RandomAccessSample] = {HoldRest};
RandomAccessSample[n_, arg_, iter__] := Block[{
iterArgs, iterLists
},
(* Separate the variables we're iterating over from the corresponding lists,
and convert range specifications into lists for later. *)

iterArgs = First /@ {iter};
iterLists = Table[
If[Head[i[[2]]] === List, i[[2]], Range @@ i[[2 ;;]]],
{i, {iter}}
];

Table[
arg /. Thread[iterArgs -> RandomChoice /@ iterLists],
n
]
];


The above can easily be parallelized by replace Table with ParallelTable and the results remain the same.

A quick check that it samples uniformly (though note here that it is repeating elements as it produces 100'000 samples from 1000 possibilities):

RandomAccessSample[100000,
x,
{x, 0, 10},
{y, Range[10]},
{z, 10}
] // Histogram