www.digitalmars.com         C & C++   DMDScript  

digitalmars.D - [OT] Algorithm question

reply "H. S. Teoh via Digitalmars-d" <digitalmars-d puremagic.com> writes:
I'm been thinking about the following problem, and thought I'd pick the
brains of the bright people around these parts...

Given a set A of n elements (let's say it's a random-access range of
size n, where n is relatively large), and a predicate P(x) that
specifies some subset of A of elements that we're interested in, what's
the best algorithm (in terms of big-O time complexity) for selecting a
random element x satisfying P(x), such that elements that satisfy P(x)
have equal probability of being chosen? (Elements that do not satisfy
P(x) are disregarded.)

Which elements of A satisfy P(x) is not known ahead of time, nor is the
relative proportion of elements that satisfy P(x) or not.

Note that A is considered to be immutable (swapping elements or
otherwise changing A is not allowed).

So far, I have come up with the following algorithms:

1) "Random shooting":

	auto r = ... /* elements of A */
	for (;;) {
		auto i = uniform(0, r.length);
		if (P(r[i])) return r[i];
	}

   Advantages: If a large percentage of elements in A satisfy P(x), then
   the loop will terminate within a small number of iterations; if the
   majority of elements satisfy P(x), this will terminate well below n
   iterations.

   Disadvantages: If only a small percentage of elements satisfy P(x),
   then this loop could take arbitrarily long to terminate, and it will
   not terminate at all if no elements satisfy P(x).

2) "Turtle walk":

	auto r = ... /* elements of A */
	int nSatisfy = 0;
	ElementType!A currentChoice;
	bool found = false;
	foreach (e; r) {
		if (P(e)) {
			if (uniform(0, nSatisfy) == 0)
			{
				currentChoice = e;
				found = true;
			}
			nSatisfy++;
		}
	}
	if (found) return currentChoice;
	else ... /* no elements satisfy P(x) */

   Advantages: If there is any element at all in A that satisfies P(x),
   it will be found. The loop always terminates after n iterations.

   Disadvantages: Even if 99% of elements in A satisfies P(x), we still
   have to traverse the entire data set before we terminate. (We cannot
   terminate before that, otherwise the probability of elements
   satisfying P(x) being selected will not be even.)

3) Permutation walk:

	auto r = ... /* elements of A */
	foreach (i; iota(0 .. r.length).randomPermutation) {
		if (P(r[i])) return r[i];
	}
	/* no elements satisfy P(x) */

   Advantages: if an element that satisfies P(x) is found early, the
   loop will terminate before n iterations. This seems like the best of
   both worlds of (1) and (2), except:

   Disadvantages: AFAIK, generating a random permutation of indices from
   0 .. n requires at least O(n) time, so any advantage we may have had
   seems to be negated.

Is there an algorithm for *incrementally* generating a random
permutation of indices? If there is, we could use that in (3) and thus
achieve the best of both worlds between early termination if an element
satisfying P(x) is found, and guaranteeing termination after n
iterations if no elements satisfying P(x) exists.


T

-- 
The early bird gets the worm. Moral: ewww...
Apr 30
next sibling parent reply Jack Stouffer <jack jackstouffer.com> writes:
On Monday, 1 May 2017 at 04:15:35 UTC, H. S. Teoh wrote:
 2) "Turtle walk":
I'd just like to point out that this algorithm doesn't satisfy
such that elements that satisfy P(x) have equal probability of 
being chosen
as the first element found in A will be chosen, and then each subsequent element has a decreasing probability of replacing the first element of P = 1/nSatisfy.
Apr 30
parent "H. S. Teoh via Digitalmars-d" <digitalmars-d puremagic.com> writes:
On Mon, May 01, 2017 at 05:03:17AM +0000, Jack Stouffer via Digitalmars-d wrote:
 On Monday, 1 May 2017 at 04:15:35 UTC, H. S. Teoh wrote:
 2) "Turtle walk":
I'd just like to point out that this algorithm doesn't satisfy
 such that elements that satisfy P(x) have equal probability of being
 chosen
as the first element found in A will be chosen, and then each subsequent element has a decreasing probability of replacing the first element of P = 1/nSatisfy.
Actually, this is exactly why it *does* satisfy the requirement. If there is exactly 1 element that satisfies P(x), then the probability of it being picked is 1. If there are 2 elements, the first element is picked with probability 1 initially, but the 2nd element has a 1/2 chance of replacing it, meaning the overall probability distribution is 1/2 and 1/2. If there are 3 elements, then by the time the 2nd element is processed, the first 2 elements are equally likely to be currently chosen (with probability 1/2 for each); when the 3rd element is found, it has 1/3 chance of replacing the previous choice, meaning there is a 2/3 chance the previous choice prevails. In that case, the 2/3 chance is equally distributed between the first two elements (they were picked with 1/2 and 1/2 probability), so the effective probability is 1/3 each after the 3rd element is processed. In general, by the time you process the n'th element, the previous elements would have had been picked with 1/(n-1) chance each, and the n'th element would be picked with 1/n chance, meaning there is a (n-1)/n chance the previous choice prevails. That (n-1)/n chance is equally distributed among the previous (n-1) elements, so effectively they are picked with 1/n chance each. The key here is that the probability that the n'th element is *not* picked is exactly (n-1)/n, which, by inductive hypothesis, must be evenly distributed among the previous (n-1) candidates, thereby making their *eventual* probability of remaining as the chosen element exactly 1/n. T -- Two wrongs don't make a right; but three rights do make a left...
Apr 30
prev sibling next sibling parent Ola Fosheim =?UTF-8?B?R3LDuHN0YWQ=?= writes:
On Monday, 1 May 2017 at 04:15:35 UTC, H. S. Teoh wrote:
 Which elements of A satisfy P(x) is not known ahead of time, 
 nor is the
 relative proportion of elements that satisfy P(x) or not.
O(N) given P(x)===false for all x... What you want is probably average case analysis, not worst case?
May 01
prev sibling next sibling parent reply John Colvin <john.loughran.colvin gmail.com> writes:
On Monday, 1 May 2017 at 04:15:35 UTC, H. S. Teoh wrote:
 I'm been thinking about the following problem, and thought I'd 
 pick the brains of the bright people around these parts...

 Given a set A of n elements (let's say it's a random-access 
 range of
 size n, where n is relatively large), and a predicate P(x) that
 specifies some subset of A of elements that we're interested 
 in, what's
 the best algorithm (in terms of big-O time complexity) for 
 selecting a
 random element x satisfying P(x), such that elements that 
 satisfy P(x)
 have equal probability of being chosen? (Elements that do not 
 satisfy
 P(x) are disregarded.)

 Which elements of A satisfy P(x) is not known ahead of time, 
 nor is the
 relative proportion of elements that satisfy P(x) or not.

 Note that A is considered to be immutable (swapping elements or
 otherwise changing A is not allowed).

 So far, I have come up with the following algorithms:

 1) "Random shooting":

 	auto r = ... /* elements of A */
 	for (;;) {
 		auto i = uniform(0, r.length);
 		if (P(r[i])) return r[i];
 	}

    Advantages: If a large percentage of elements in A satisfy 
 P(x), then
    the loop will terminate within a small number of iterations; 
 if the
    majority of elements satisfy P(x), this will terminate well 
 below n
    iterations.

    Disadvantages: If only a small percentage of elements 
 satisfy P(x),
    then this loop could take arbitrarily long to terminate, and 
 it will
    not terminate at all if no elements satisfy P(x).
You can recover O(n) calls to P by caching the results. You can recover O(n) for both P and rng by progressively removing elements from r or any other method that results in an easily samplable set of all elements of r not yet visited (and therefore possible candidates). I like this one (warning, just thought of, untested): auto r = ... /* elements of A */ auto nRemaining = r.length; while (nRemaining) { auto i = uniform(0, nRemaining); if (P(r[i])) return r[i]; else swap(r[i], r[--nRemaining]); }
 3) Permutation walk:

 	auto r = ... /* elements of A */
 	foreach (i; iota(0 .. r.length).randomPermutation) {
 		if (P(r[i])) return r[i];
 	}
 	/* no elements satisfy P(x) */

    Advantages: if an element that satisfies P(x) is found 
 early, the
    loop will terminate before n iterations. This seems like the 
 best of
    both worlds of (1) and (2), except:

    Disadvantages: AFAIK, generating a random permutation of 
 indices from
    0 .. n requires at least O(n) time, so any advantage we may 
 have had
    seems to be negated.

 Is there an algorithm for *incrementally* generating a random 
 permutation of indices? If there is, we could use that in (3) 
 and thus achieve the best of both worlds between early 
 termination if an element satisfying P(x) is found, and 
 guaranteeing termination after n iterations if no elements 
 satisfying P(x) exists.


 T
Yes. As a matter of fact it would be the logical extension of my algorithm above for when you can't or don't want to change r (again, only just thought of this, untested...): auto r = ... /* elements of A */ auto indices = iota(r.length).array; auto nRemaining = r.length; while (nRemaining) { auto i = uniform(0, nRemaining); if (P(r[indices[i]])) return r[indices[i]]; else swap(indices[i], indices[--nRemaining]); } You could also do the same thing but with an array of pointers to elements of r instead of indices.
May 01
next sibling parent reply Ola Fosheim =?UTF-8?B?R3LDuHN0YWQ=?= writes:
On Monday, 1 May 2017 at 08:42:05 UTC, John Colvin wrote:
 I like this one (warning, just thought of, untested):

 auto r = ... /* elements of A */
 auto nRemaining = r.length;
 while (nRemaining) {
 	auto i = uniform(0, nRemaining);
 	if (P(r[i])) return r[i];
 	else swap(r[i], r[--nRemaining]);
 }
Yes, this is the standard text-book way of doing it, still O(N) of course. The other standard-text-book way is to generate an array of indexes and mutate that instead, still O(N). If you cache in a heap you get O(N log N). Anyway, this kind of worst case analysis doesn't really help. And neither does average case, which will be O(1) assuming 50% match P. So you need is to specify what kind of average case analysis you want. For instance expressed as f(N,S) where N is the number of elements and S is the number of elements that satisfies P.
May 01
parent Ola Fosheim =?UTF-8?B?R3LDuHN0YWQ=?= writes:
On Monday, 1 May 2017 at 09:01:33 UTC, Ola Fosheim Grøstad wrote:
 array of indexes and mutate that instead, still O(N). If you 
 cache in a heap you get O(N log N).
To avoid confusion: I didn't mean a literal "heap", but some kind of search structure that tracks and draws from unused indice ranges in a way that does not require O(N) initialization and can be implemented as a single array on the stack. I don't remember what it is called, but you draw on each node based on the number of children (kind of like the ones used in markov-chain type situations). I assume this can be maintained in O(log N) per iteration with O(1) initialization. E.g. 1. push_indices(0,N-1) //data (0,N-1) 2. i = draw_and_remove_random_index() // data could now be (0,i-1), (i+1,N-1) 3. if( i == -1 ) return 4. dostuff(i) 5. goto 2 Even though it is O(N log N) it could outperform other solutions if you only want to draw a small number of elements. And a small array with a linear search (tracking used indices) would outperform that if you only want to draw say 8 elements. If you want to draw close to N elements then you would be better off just randomizing a full array (with elements from 0 to N-1). Or you could do all of these, start with linear search, then convert to a tree like search, then convert to full enumeration. In general I don't think asymptotical analysis will do much good. I think you need to look on actual costs for typical N, P and different distributions of satisfied P, e.g. cost of loop iteration: 10 average cost of P(x)=>false: 1 average cost of P(x)=>true: 1000 probability for P(a[i]) for each i
May 01
prev sibling parent "H. S. Teoh via Digitalmars-d" <digitalmars-d puremagic.com> writes:
On Mon, May 01, 2017 at 08:42:05AM +0000, John Colvin via Digitalmars-d wrote:
[...]
 You can recover O(n) calls to P by caching the results.
Sorry, I forgot to specify that P(x) can be different each time. Also, the input A may also be different each time, albeit remain the same length (though the expectation is that it will be mostly the same -- not that it matters, though, I don't think any of the proposed solutions here count on that).
 You can recover O(n) for both P and rng by progressively removing
 elements from r or any other method that results in an easily
 samplable set of all elements of r not yet visited (and therefore
 possible candidates).
[...] Yes, I realize that if r was mutable, there are existing simple solutions. The problem is that it's not mutable, or rather, should not be mutated by this selection algorithm. [...]
 Yes. As a matter of fact it would be the logical extension of my
 algorithm above for when you can't or don't want to change r (again,
 only just thought of this, untested...):
 
 auto r = ... /* elements of A */
 auto indices = iota(r.length).array;
 auto nRemaining = r.length;
 while (nRemaining) {
 	auto i = uniform(0, nRemaining);
 	if (P(r[indices[i]])) return r[indices[i]];
 	else swap(indices[i], indices[--nRemaining]);
 }
 
 You could also do the same thing but with an array of pointers to
 elements of r instead of indices.
The trouble with this is that you have to allocate an array of indices, and r.length is rather large, so creating this array is O(n) each time. But I think Ivan Kazmenko has an excellent idea to cache this index array. Since, by symmetry, it doesn't matter if the starting indices array was the identity permutation, we could just initialize this array once, and the next time just skip to the loop directly (reusing the previous (possibly partially) permuted indices as the starting point). It does require O(n) memory, which is unfortunate, but at least we only pay for that once. I was hoping for something more along the lines of a random permutation algorithm with sublinear (if not O(1)) space complexity that can return elements of the permutation consecutively on-demand, like a lazy range. Something along the lines of: auto seed = ... /* get random number */ auto lazyRange = generatePermutation(seed, n); where generatePermutation is some incremental algorithm that produces a unique permutation of 0 .. n per seed value. Ideally, this algorithm would only do as much work as is necessary to produce the first k elements of the permutation, so that if the main loop terminates early (we found an element satisfying P(x)) we don't waste time generating the rest of the permutation. T -- "No, John. I want formats that are actually useful, rather than over-featured megaliths that address all questions by piling on ridiculous internal links in forms which are hideously over-complex." -- Simon St. Laurent on xml-dev
May 01
prev sibling next sibling parent reply Ola Fosheim =?UTF-8?B?R3LDuHN0YWQ=?= writes:
On Monday, 1 May 2017 at 04:15:35 UTC, H. S. Teoh wrote:
 Is there an algorithm for *incrementally* generating a random 
 permutation of indices?
Does it exist? Yes, because you can build permutation tables for it, but you'll have to find the closed form for it that is efficient... Can you do that without selecting a specific N? It is easy for N=2 and N=3 at least. E.g. For array length 2 you get the 2 permutations: 0 1, 1 0 Then you just select one of them at random (you know the number of permutations) So if it exists you'll should probably do a search for permutations. "How do I efficiently generate permutation number N" Of course for smaller N you could generate a big array of bytes and compress it by having arrays of slices into it (as many permutations will share index sequences).
May 01
parent reply Ola Fosheim =?UTF-8?B?R3LDuHN0YWQ=?= writes:
On Monday, 1 May 2017 at 08:44:25 UTC, Ola Fosheim Grøstad wrote:
 Does it exist? Yes, because you can build permutation tables 
 for it, but you'll have to find the closed form for it that is 
 efficient... Can you do that without selecting a specific N? It 
 is easy for N=2 and N=3 at least.

 E.g.

 For array length 2 you get the 2 permutations: 0 1, 1 0
 Then you just select one of them at random (you know the number 
 of permutations)

 So if it exists you'll should probably do a search for 
 permutations. "How do I efficiently generate permutation number 
 N"

 Of course for smaller N you could generate a big array of bytes 
 and compress it by having arrays of slices into it (as many 
 permutations will share index sequences).
Keep also in mind that you can split the original array in two, so it might be sufficient to have the permutations for various 2^M sizes if those are easier to generate (my guess, maybe through some xor magic?) E.g. if N = 21 then it can be split into 16 + (4 + 1) So you can draw like this: 1. randomly select group_16 over group_5 based with 16/21 probability) 2. if group_5 selected: ramdomly select group_4 over group 1 on 4/5 probability etc. That way you only need log(N) permutation sequences to keep track of. Which for all practical purposes are going to stay pretty low (<20?) Right?
May 01
parent reply Ola Fosheim =?UTF-8?B?R3LDuHN0YWQ=?= writes:
On Monday, 1 May 2017 at 10:51:46 UTC, Ola Fosheim Grøstad wrote:
 Keep also in mind that you can split the original array in two, 
 so it might be sufficient to  have the permutations for various 
 2^M sizes if those are easier to generate (my guess, maybe 
 through some xor magic?)

 E.g. if N = 21 then it can be split into 16 + (4 + 1)

 So you can draw like this:

 1. randomly select group_16 over group_5 based with 16/21 
 probability)

 2. if group_5 selected: ramdomly select group_4 over group 1 on 
 4/5 probability

 etc.

 That way you only need log(N) permutation sequences to keep 
 track of. Which for all practical purposes are going to stay 
 pretty low (<20?)

 Right?
But if you don't use a hardware random generator and are happy with more sloppy pseudo-randomness then you would be better of ditching the whole concept of permutations and replace it with some cylic random generators intead using the same technique I just outlined. E.g. find a set of cyclic random generators and break down N into a sum of these cycle-sizes, then just keep track of the seed for each. If they are 2^N then you could use xor to get more randomness between runs. Also in the algorithms above you need to change the probabilities each time you take away one index from a group (you don't want to draw from an empty group).
May 01
parent reply Ola Fosheim =?UTF-8?B?R3LDuHN0YWQ=?= writes:
On Monday, 1 May 2017 at 11:08:56 UTC, Ola Fosheim Grøstad wrote:
 E.g. find a set of cyclic random generators and break down N 
 into a sum of these cycle-sizes, then just keep track of the 
 seed for each. If they are 2^N then you could use xor to get 
 more randomness between runs.

 Also in the algorithms above you need to change the 
 probabilities each time you take away one index from a group 
 (you don't want to draw from an empty group).
Well, actually, just select the single cyclic generator that is larger or equal to N, then just redraw if it returns a value >= N. Duh! Sloppy thinking. I would think you should be able to find some prime sized ones that will get the next index with an insignificant number of redraws. But permutations is the way to go if you want scientific quality.
May 01
parent reply Ola Fosheim =?UTF-8?B?R3LDuHN0YWQ=?= writes:
On Monday, 1 May 2017 at 11:15:28 UTC, Ola Fosheim Grøstad wrote:
 But permutations is the way to go if you want scientific 
 quality.
I take that back: «Polynomial pseudo-random number generator via cyclic phase» http://www.sciencedirect.com/science/article/pii/S0378475409001463 Might be what you are looking for and if not it probably references other papers that fits the bill.
May 01
parent reply Ola Fosheim =?UTF-8?B?R3LDuHN0YWQ=?= writes:
Ah, found a sensible search term: «pseudorandom permutation»

http://cse.iitkgp.ac.in/~debdeep/courses_iitkgp/FCrypto/slides/LubyRackoff.pdf

If you go down this path, please share links. Useful stuff 
actually.
May 01
parent Ola Fosheim =?UTF-8?B?R3LDuHN0YWQ=?= writes:
On Monday, 1 May 2017 at 11:57:32 UTC, Ola Fosheim Grøstad wrote:
 Ah, found a sensible search term: «pseudorandom permutation»

 http://cse.iitkgp.ac.in/~debdeep/courses_iitkgp/FCrypto/slides/LubyRackoff.pdf

 If you go down this path, please share links. Useful stuff 
 actually.
A readable paper for most people: http://www.pcg-random.org/pdf/toms-oneill-pcg-family-v1.02.pdf «some users would prefer a much stronger property analogous to uniformity: that over the full period of the generator, every possible k-tuple will occur, and it will occur the same number of times.» PCG claims to have arbitrary cycle sizes. I haven't looked at it, but it is available as C/C++.
May 01
prev sibling next sibling parent reply Moritz Maxeiner <moritz ucworks.org> writes:
On Monday, 1 May 2017 at 04:15:35 UTC, H. S. Teoh wrote:
 Given a set A of n elements (let's say it's a random-access 
 range of
 size n, where n is relatively large), and a predicate P(x) that
 specifies some subset of A of elements that we're interested 
 in, what's
 the best algorithm (in terms of big-O time complexity) for 
 selecting a
 random element x satisfying P(x), such that elements that 
 satisfy P(x)
 have equal probability of being chosen? (Elements that do not 
 satisfy
 P(x) are disregarded.)

 Which elements of A satisfy P(x) is not known ahead of time, 
 nor is the
 relative proportion of elements that satisfy P(x) or not.
Works only with a random access range and I haven't done the full analysis (and I'm a bit rusty on probability), but my first thought was random divide and conquer: ElementType!A* select(A,P)(A r) { // Recursion anchor if (r.length == 1) { if (P(r[0])) return r[0]; else return null; // Recurse randomly with p=0.5 into either the left, or right half of r } else { ElementType!A* e; ElementType!A[][2] half; half[0] = r[0..floor(r.length/2)]; half[1] = r[ceil(r.length/2)..$]; ubyte coinFlip = uniform(0,1) > 0; // Recurse in one half and terminate if e is found there e = select(half[coinFlip]); if (e) return e; // Otherwise, recurse into other half return select(half[1 - coinFlip]); } } As stated above, I haven't done the full analysis, but intuitively speaking (which can be wrong, of course), the p=0.5 recursion ought satisfy the condition of elements satisfying P(x) being chosen uniformly; also intuitively, I'd expect the expected runtime for a uniform distribution of elements satisfying P(x) to be around O(log N). Worst-case would be that it has to inspect every element in r once => O(N)
May 01
parent reply Timon Gehr <timon.gehr gmx.ch> writes:
On 01.05.2017 11:51, Moritz Maxeiner wrote:
 On Monday, 1 May 2017 at 04:15:35 UTC, H. S. Teoh wrote:
 Given a set A of n elements (let's say it's a random-access range of
 size n, where n is relatively large), and a predicate P(x) that
 specifies some subset of A of elements that we're interested in, what's
 the best algorithm (in terms of big-O time complexity) for selecting a
 random element x satisfying P(x), such that elements that satisfy P(x)
 have equal probability of being chosen? (Elements that do not satisfy
 P(x) are disregarded.)

 Which elements of A satisfy P(x) is not known ahead of time, nor is the
 relative proportion of elements that satisfy P(x) or not.
Works only with a random access range and I haven't done the full analysis (and I'm a bit rusty on probability), but my first thought was random divide and conquer: ElementType!A* select(A,P)(A r) { // Recursion anchor if (r.length == 1) { if (P(r[0])) return r[0]; else return null; // Recurse randomly with p=0.5 into either the left, or right half of r } else { ElementType!A* e; ElementType!A[][2] half; half[0] = r[0..floor(r.length/2)]; half[1] = r[ceil(r.length/2)..$]; ubyte coinFlip = uniform(0,1) > 0;
This deterministically chooses 0. (uniform is right-exclusive.) I assume you mean uniform(0,2).
     // Recurse in one half and terminate if e is found there
     e = select(half[coinFlip]);
     if (e) return e;
     // Otherwise, recurse into other half
     return select(half[1 - coinFlip]);
   }
 }

 As stated above, I haven't done the full analysis, but intuitively
 speaking (which can be wrong, of course), the p=0.5 recursion ought
 satisfy the condition of elements satisfying P(x) being chosen
 uniformly;
Counterexample: [1,0,1,1]. The first element is chosen with probability 1/2, which is not 1/3.
 also intuitively, I'd expect the expected runtime for a
 uniform distribution of elements satisfying P(x)  to be around O(log N).
I don't understand this input model.
 Worst-case would be that it has to inspect every element in r once => O(N)
May 01
parent Moritz Maxeiner <moritz ucworks.org> writes:
On Monday, 1 May 2017 at 09:58:39 UTC, Timon Gehr wrote:
 On 01.05.2017 11:51, Moritz Maxeiner wrote:
     [...]
This deterministically chooses 0. (uniform is right-exclusive.) I assume you mean uniform(0,2).
Yes.
 [...]
Counterexample: [1,0,1,1]. The first element is chosen with probability 1/2, which is not 1/3.
Dang. Serves me right for not doing the full analysis. Thanks for the counter example and sorry for wasting your time.
 [...]
I don't understand this input model.
Since the idea was nonsense, it doesn't matter, anyway.
May 01
prev sibling next sibling parent reply Mike B Johnson <Mikey Ikes.com> writes:
On Monday, 1 May 2017 at 04:15:35 UTC, H. S. Teoh wrote:
 I'm been thinking about the following problem, and thought I'd 
 pick the brains of the bright people around these parts...

 [...]
Since most real world problems would require selecting elements more than once it may be far more efficient to sort by P(x)(filter) then simply select a random element. e.g., a = A.filter(x->P(x)) // Creates a new set a where only the elements of A that satisfy P(x) are added ... e = a.random; this is O(1) if you only have to filter once(you can create a container that always "sorts" on P(x) so to speak.. like a sorted dictionary).
May 01
parent "H. S. Teoh via Digitalmars-d" <digitalmars-d puremagic.com> writes:
On Mon, May 01, 2017 at 11:32:19AM +0000, Mike B Johnson via Digitalmars-d
wrote:
[...]
 Since most real world problems would require selecting elements more
 than once it may be far more efficient to sort by P(x)(filter) then
 simply select a random element.
 
 e.g.,
 
 a = A.filter(x->P(x))  // Creates a new set a where only the elements
 of A that satisfy P(x) are added
 
 ...
 
 e = a.random;
 
 
 this is O(1) if you only have to filter once(you can create a
 container that always "sorts" on P(x) so to speak.. like a sorted
 dictionary).
[...] Sorry, I forgot to specify that P(x) may change each time, as may the input A. So caching would be of little help in this case. T -- This is a tpyo.
May 01
prev sibling next sibling parent reply Andrea Fontana <nospam example.com> writes:
On Monday, 1 May 2017 at 04:15:35 UTC, H. S. Teoh wrote:
 1) "Random shooting":

 	auto r = ... /* elements of A */
 	for (;;) {
 		auto i = uniform(0, r.length);
 		if (P(r[i])) return r[i];
 	}

    Advantages: If a large percentage of elements in A satisfy 
 P(x), then
    the loop will terminate within a small number of iterations; 
 if the
    majority of elements satisfy P(x), this will terminate well 
 below n
    iterations.

    Disadvantages: If only a small percentage of elements 
 satisfy P(x),
    then this loop could take arbitrarily long to terminate, and 
 it will
    not terminate at all if no elements satisfy P(x).
If you find an element that doesn't satisfy P(x) move it on top of array (swap!) and then try again with uniform(1, r.length - 1) and so on. It terminates in this way. Andrea
May 01
next sibling parent reply Andrea Fontana <nospam example.com> writes:
On Monday, 1 May 2017 at 12:31:48 UTC, Andrea Fontana wrote:
 If you find an element that doesn't satisfy P(x) move it on top 
 of array (swap!) and then try again with uniform(1, r.length - 
 1) and so on.
Whoops i mean uniform(1, r.length) of course :)
May 01
parent Ola Fosheim =?UTF-8?B?R3LDuHN0YWQ=?= writes:
On Monday, 1 May 2017 at 12:38:32 UTC, Andrea Fontana wrote:
 On Monday, 1 May 2017 at 12:31:48 UTC, Andrea Fontana wrote:
 If you find an element that doesn't satisfy P(x) move it on 
 top of array (swap!) and then try again with uniform(1, 
 r.length - 1) and so on.
Whoops i mean uniform(1, r.length) of course :)
I guess you meant to the end, but if mutating the array is ok, then you don't have to swap. Just move in the last element and chop 1 off the length. Another cache-friendly mutation solution (assuming PRNG is cheap) is to have two phases: Phase 1: linear walkover the array and select elements with P(a[i])==true with probability 1/N the overwrite the ones with P(a[i]) = false. when done truncate the array. Phase 2: regular random selection from the array were all elements satisfy P. The OP underspecified the requirements... There is a gazillion variations... :-P
May 01
prev sibling parent "H. S. Teoh via Digitalmars-d" <digitalmars-d puremagic.com> writes:
On Mon, May 01, 2017 at 12:31:48PM +0000, Andrea Fontana via Digitalmars-d
wrote:
 On Monday, 1 May 2017 at 04:15:35 UTC, H. S. Teoh wrote:
 1) "Random shooting":
 
 	auto r = ... /* elements of A */
 	for (;;) {
 		auto i = uniform(0, r.length);
 		if (P(r[i])) return r[i];
 	}
 
    Advantages: If a large percentage of elements in A satisfy P(x), then
    the loop will terminate within a small number of iterations; if the
    majority of elements satisfy P(x), this will terminate well below n
    iterations.
 
    Disadvantages: If only a small percentage of elements satisfy P(x),
    then this loop could take arbitrarily long to terminate, and it will
    not terminate at all if no elements satisfy P(x).
 
If you find an element that doesn't satisfy P(x) move it on top of array (swap!) and then try again with uniform(1, r.length - 1) and so on. It terminates in this way.
[...] The problem is that swapping elements is not allowed for this particular problem. If swapping was allowed this would be a much easier problem to solve. :-) T -- Chance favours the prepared mind. -- Louis Pasteur
May 01
prev sibling next sibling parent reply Ivan Kazmenko <gassa mail.ru> writes:
On Monday, 1 May 2017 at 04:15:35 UTC, H. S. Teoh wrote:
 Given a set A of n elements (let's say it's a random-access 
 range of
 size n, where n is relatively large), and a predicate P(x) that
 specifies some subset of A of elements that we're interested 
 in, what's
 the best algorithm (in terms of big-O time complexity) for 
 selecting a
 random element x satisfying P(x), such that elements that 
 satisfy P(x)
 have equal probability of being chosen? (Elements that do not 
 satisfy
 P(x) are disregarded.)
I'd like to note here that, if you make use of the same P(x) many times (instead of different predicates on each call), it makes sense to spend O(n) time and memory filtering by that predicate and storing the result, and then answer each query in O(1).
 3) Permutation walk:

 	auto r = ... /* elements of A */
 	foreach (i; iota(0 .. r.length).randomPermutation) {
 		if (P(r[i])) return r[i];
 	}
 	/* no elements satisfy P(x) */

    Advantages: if an element that satisfies P(x) is found 
 early, the
    loop will terminate before n iterations. This seems like the 
 best of
    both worlds of (1) and (2), except:

    Disadvantages: AFAIK, generating a random permutation of 
 indices from
    0 .. n requires at least O(n) time, so any advantage we may 
 have had
    seems to be negated.

 Is there an algorithm for *incrementally* generating a random 
 permutation of indices? If there is, we could use that in (3) 
 and thus achieve the best of both worlds between early 
 termination if an element satisfying P(x) is found, and 
 guaranteeing termination after n iterations if no elements 
 satisfying P(x) exists.
Yes, there is. There are actually two variations of Fisher-Yates shuffle: (https://en.wikipedia.org/wiki/Fisher%E2%80%93Yates_shuffle) 1. auto p = n.iota.array; foreach (pos; 0..n) { auto otherPos = uniform (0, pos + 1); swap (p[pos], p[otherPos]); } When we look at this after k-th iteration, the first k elements are randomly and uniformly permuted, and the rest (n-k) are left untouched. 2. auto p = n.iota.array; foreach (pos; 0..n) { auto otherPos = uniform (pos, n); swap (p[pos], p[otherPos]); } When we look at this after k-th iteration, the first k elements are a random combination of all n elements, and this combination is randomly and uniformly permuted. So, the second variation is what we need: each new element is randomly and uniformly selected from all the elements left. Once we get the first element satisfying the predicate, we can just terminate the loop. If there are m out of n elements satisfying the predicate, the average number of steps is n/m. Now, the problem is that both of these allocate n "size_t"-s of memory to start with. And your problem does not allow to shuffle the elements of the original array in place, so we do need an external permutation for these algorithms. However, there are at least two ways to mitigate that: (I) We can allocate the permutation once using n time and memory, and then, on every call, just reuse it in its current state in n/m time. It does not matter if the permutation is not identity permutation: by symmetry argument, any starting permutation will do just fine. (II) We can store the permutation p in an associative array instead of a regular array, actually storing only the elements accessed at least once, and assuming other elements to satisfy the identity p[x] = x. So, if we finish in n/m steps on average, the time and extra memory used will be O(n/m) too. I can put together an example implementation if this best satisfies your requirements. Ivan Kazmenko.
May 01
parent "H. S. Teoh via Digitalmars-d" <digitalmars-d puremagic.com> writes:
On Mon, May 01, 2017 at 02:38:09PM +0000, Ivan Kazmenko via Digitalmars-d wrote:
 On Monday, 1 May 2017 at 04:15:35 UTC, H. S. Teoh wrote:
 Given a set A of n elements (let's say it's a random-access range of
 size n, where n is relatively large), and a predicate P(x) that
 specifies some subset of A of elements that we're interested in,
 what's the best algorithm (in terms of big-O time complexity) for
 selecting a random element x satisfying P(x), such that elements
 that satisfy P(x) have equal probability of being chosen? (Elements
 that do not satisfy P(x) are disregarded.)
I'd like to note here that, if you make use of the same P(x) many times (instead of different predicates on each call), it makes sense to spend O(n) time and memory filtering by that predicate and storing the result, and then answer each query in O(1).
Unfortunately, P(x) could be different each time (sorry, should have stated this explicitly), and A may also change in between calls to this random selection algorithm. So filtering would not be a good approach. [...]
 There are actually two variations of Fisher-Yates shuffle:
 (https://en.wikipedia.org/wiki/Fisher%E2%80%93Yates_shuffle)
 
 1.
     auto p = n.iota.array;
     foreach (pos; 0..n) {
         auto otherPos = uniform (0, pos + 1);
         swap (p[pos], p[otherPos]);
     }
 
 When we look at this after k-th iteration, the first k elements are
 randomly and uniformly permuted, and the rest (n-k) are left
 untouched.
 
 2.
     auto p = n.iota.array;
     foreach (pos; 0..n) {
         auto otherPos = uniform (pos, n);
         swap (p[pos], p[otherPos]);
     }
 
 When we look at this after k-th iteration, the first k elements are a
 random combination of all n elements, and this combination is randomly
 and uniformly permuted.  So, the second variation is what we need:
 each new element is randomly and uniformly selected from all the
 elements left.  Once we get the first element satisfying the
 predicate, we can just terminate the loop.  If there are m out of n
 elements satisfying the predicate, the average number of steps is n/m.
But wouldn't creating the array p already be O(n)? Albeit, somewhat faster than n iterations of the main loop since copying iota ought to be somewhat cheaper than generating a random number and swapping two elements.
 Now, the problem is that both of these allocate n "size_t"-s of memory
 to start with.  And your problem does not allow to shuffle the
 elements of the original array in place, so we do need an external
 permutation for these algorithms.  However, there are at least two
 ways to mitigate that:
 
 (I)
 We can allocate the permutation once using n time and memory, and
 then, on every call, just reuse it in its current state in n/m time.
 It does not matter if the permutation is not identity permutation: by
 symmetry argument, any starting permutation will do just fine.
I like this idea very much. This lets us only pay once for creating the index array, and reuse it almost "for free" thereafter. Still, the O(n) space requirement could potentially be onerous, since n is expected to be rather large.
 (II)
 We can store the permutation p in an associative array instead of a
 regular array, actually storing only the elements accessed at least
 once, and assuming other elements to satisfy the identity p[x] = x.
 So, if we finish in n/m steps on average, the time and extra memory
 used will be O(n/m) too.  I can put together an example implementation
 if this best satisfies your requirements.
[...] This is interesting, and relaxes the O(n) space requirement initially. After enough queries, though, I'd expect the AA to eventually grow to n elements (as there is a chance some queries will end up finding no elements that satisfy P(x) so it would generate the entire permutation). But I like your idea of reusing the index array. It's definitely a possibility! T -- Your inconsistency is the only consistent thing about you! -- KD
May 01
prev sibling parent reply MysticZach <reachzach ggmail.com> writes:
On Monday, 1 May 2017 at 04:15:35 UTC, H. S. Teoh wrote:
 Given a set A of n elements (let's say it's a random-access 
 range of
 size n, where n is relatively large), and a predicate P(x) that
 specifies some subset of A of elements that we're interested 
 in, what's
 the best algorithm (in terms of big-O time complexity) for 
 selecting a
 random element x satisfying P(x), such that elements that 
 satisfy P(x)
 have equal probability of being chosen? (Elements that do not 
 satisfy
 P(x) are disregarded.)
Here's how I would do it: // returns a random index of array satisfying P(x), -1 if not found int randomlySatisfy(A[] array) { if (array.length == 0) return -1; int i = uniform(0, array.length); return randomlySatisfyImpl(array, i); } // recursive function private int randomlySatisfyImpl(A[] da, int i) { if (P(da[i])) return i; if (da.length == 1) return -1; // choose a random partition proportionally auto j = uniform(da.length - 1); if (j < i) { // the lower partition int a = randomlySatisfyImpl(da[0..i], j); if (a != -1) return a; else return randomlySatisfyImpl(da[i+1 .. da.length], j - (i + 1)); } else { // higher partition, investigate in reverse order int a = randomlySatisfyImpl(da[i+1 .. da.length], j - (i + 1)); if (a != -1) return i +1 + a; else return i + 1 + randomlySatisfyImpl(da[0..i], j); } } The goal is to have the first hit be the one you return. The method: if a random pick doesn't satisfy, randomly choose the partition greater than or less than based on uniform(0..array.length-1), and do the same procedure on that partition, reusing the random index to avoid having to call uniform twice on each recursion (one to choose a partition and one to choose an index within that partition). If the probability of choosing a partition is precisely proportional to the number of elements in that partition, it seems to me that the result will be truly random, but I'm not sure.
May 01
next sibling parent MysticZach <reachzach ggmail.com> writes:
On Monday, 1 May 2017 at 16:56:58 UTC, MysticZach wrote:
    // choose a random partition proportionally
    auto j = uniform(da.length - 1);
    if (j < i) {
       // the lower partition
       int a = randomlySatisfyImpl(da[0..i], j);
       if (a != -1) return a;
       else return randomlySatisfyImpl(da[i+1 .. da.length], j - 
 (i + 1));
    }
    else {
       // higher partition, investigate in reverse order
       int a = randomlySatisfyImpl(da[i+1 .. da.length], j - (i 
 + 1));
       if (a != -1) return i +1 + a;
       else return i + 1 + randomlySatisfyImpl(da[0..i], j);
The line above has a bug. Replace it with: else { a = randomlySatisfyImpl(da[0..i], j); return (a == -1) ? -1 : i + 1 + a; }
    }
 }
But the idea's the same. Hopefully it's clear.
May 01
prev sibling parent reply MysticZach <reachzach ggmail.com> writes:
On Monday, 1 May 2017 at 16:56:58 UTC, MysticZach wrote:
 The goal is to have the first hit be the one you return. The 
 method: if a random pick doesn't satisfy, randomly choose the 
 partition greater than or less than based on 
 uniform(0..array.length-1), and do the same procedure on that 
 partition, reusing the random index to avoid having to call 
 uniform twice on each recursion (one to choose a partition and 
 one to choose an index within that partition). If the 
 probability of choosing a partition is precisely proportional 
 to the number of elements in that partition, it seems to me 
 that the result will be truly random, but I'm not sure.
Now I'm questioning this, because if the positive cases are unevenly distributed, i.e., [111111000100], the last one has about 50% chance to get picked instead of a 1 in 7 chance with my method. I guess you'd need to build a whole new array like the others are saying.
May 01
parent reply Ivan Kazmenko <gassa mail.ru> writes:
On Monday, 1 May 2017 at 21:54:43 UTC, MysticZach wrote:
 On Monday, 1 May 2017 at 16:56:58 UTC, MysticZach wrote:
 The goal is to have the first hit be the one you return. The 
 method: if a random pick doesn't satisfy, randomly choose the 
 partition greater than or less than based on 
 uniform(0..array.length-1), and do the same procedure on that 
 partition, reusing the random index to avoid having to call 
 uniform twice on each recursion (one to choose a partition and 
 one to choose an index within that partition). If the 
 probability of choosing a partition is precisely proportional 
 to the number of elements in that partition, it seems to me 
 that the result will be truly random, but I'm not sure.
Now I'm questioning this, because if the positive cases are unevenly distributed, i.e., [111111000100], the last one has about 50% chance to get picked instead of a 1 in 7 chance with my method. I guess you'd need to build a whole new array like the others are saying.
A pity; it sounded nice! But yeah, once we pick the right ~half, we have to completely traverse it before paying any attention to the left ~half. I hope some part of the idea is still salvageable. For example, what if we put the intervals in a queue instead of a stack? Ivan Kazmenko.
May 02
parent reply Ivan Kazmenko <gassa mail.ru> writes:
On Tuesday, 2 May 2017 at 10:35:46 UTC, Ivan Kazmenko wrote:
 I hope some part of the idea is still salvageable.
 For example, what if we put the intervals in a queue instead of 
 a stack?
I tried to implement a similar approach, but instead of a queue or a stack, I used a random-access array of intervals. Sadly, it is still far from uniform, since it does not take interval lengths into account, and I don't see how to do that without at least O(log n) time per interval insertion or deletion. Implementation and empiric frequencies for n=5 elements in a permutation: http://ideone.com/3zSxLN Ivan Kazmenko.
May 02
parent reply MysticZach <reachzach ggmail.com> writes:
On Tuesday, 2 May 2017 at 11:27:17 UTC, Ivan Kazmenko wrote:
 On Tuesday, 2 May 2017 at 10:35:46 UTC, Ivan Kazmenko wrote:
 I hope some part of the idea is still salvageable.
 For example, what if we put the intervals in a queue instead 
 of a stack?
I tried to implement a similar approach, but instead of a queue or a stack, I used a random-access array of intervals. Sadly, it is still far from uniform, since it does not take interval lengths into account, and I don't see how to do that without at least O(log n) time per interval insertion or deletion. Implementation and empiric frequencies for n=5 elements in a permutation: http://ideone.com/3zSxLN Ivan Kazmenko.
Well, I thought of another idea that may not be technically random, but which I imagine would get pretty close in real world use cases: 1. Test a random element in the array. If it satisfies P, return it. 2. Choose a hopper interval, h, that is relatively prime to the number of elements in the array, n. You could do this by randomly selecting from a pre-made list of primes of various orders of magnitude larger than n, with h = the prime % n. 3. Hop along the array, testing each element as you go. Increment a counter. If you reach the end of the array, cycle back to the beginning starting with the remainder of h that you didn't use. I think that h being relatively prime means it will thereby hit every element in the array. 4. Return the first hit. If the counter reaches n, return failure. The way I see it, with a random start position, and a *mostly* random interval, the only way to defeat the randomness of the result would be if the elements that satisfy P were dispersed precisely according to the random interval specified for that particular evocation of the function — in which case the first in the dispersed set would be chosen more often. This would require a rare convergence depending on both n and h, and would not repeat from one call to the next.
May 02
parent reply MysticZach <reachzach ggmail.com> writes:
On Tuesday, 2 May 2017 at 13:44:03 UTC, MysticZach wrote:
 1. Test a random element in the array. If it satisfies P, 
 return it.
 2. Choose a hopper interval, h, that is relatively prime to the 
 number of elements in the array, n. You could do this by 
 randomly selecting from a pre-made list of primes of various 
 orders of magnitude larger than n, with h = the prime % n.
 3. Hop along the array, testing each element as you go. 
 Increment a counter. If you reach the end of the array, cycle 
 back to the beginning starting with the remainder of h that you 
 didn't use. I think that h being relatively prime means it will 
 thereby hit every element in the array.
 4. Return the first hit. If the counter reaches n, return 
 failure.
Taking this a step further, it occurred to me that you could use *any* hopping interval from 1 to array.length, if the length of the array were prime. So artificially extend the array and just skip a jump when you land in the extended area. And since you'd skip lots of jumps if the extension were any bigger that the hopping interval, reduce its length to modulo the hopping interval. // returns a random index of array satisfying P(x), -1 if not found int randomlySatisfy(A[] array) { auto length = array.length; if (!length) return -1; auto i = uniform(0, length); auto hop = uniform(1, length); // greatest prime < 2^^31, for simplicity enum prime = 2147483477; assert (length <= prime); auto skipLength = ((prime - length) % hop) + length; auto count = 0; for (;;) { if (P(a[i])) return i; ++count; if (count == length) return -1; i += hop; if (i < length) continue; if (i < skipLength) i += hop; i -= skipLength; } return -1; } This solution is stupidly simple and I haven't tested it. But I believe it's truly random, since the hopping interval is arbitrary. Who knows, it might work.
May 02
parent reply Timon Gehr <timon.gehr gmx.ch> writes:
On 02.05.2017 22:42, MysticZach wrote:
 On Tuesday, 2 May 2017 at 13:44:03 UTC, MysticZach wrote:
 1. Test a random element in the array. If it satisfies P, return it.
 2. Choose a hopper interval, h, that is relatively prime to the number
 of elements in the array, n. You could do this by randomly selecting
 from a pre-made list of primes of various orders of magnitude larger
 than n, with h = the prime % n.
 3. Hop along the array, testing each element as you go. Increment a
 counter. If you reach the end of the array, cycle back to the
 beginning starting with the remainder of h that you didn't use. I
 think that h being relatively prime means it will thereby hit every
 element in the array.
 4. Return the first hit. If the counter reaches n, return failure.
Taking this a step further, it occurred to me that you could use *any* hopping interval from 1 to array.length, if the length of the array were prime. So artificially extend the array and just skip a jump when you land in the extended area. And since you'd skip lots of jumps if the extension were any bigger that the hopping interval, reduce its length to modulo the hopping interval. // returns a random index of array satisfying P(x), -1 if not found int randomlySatisfy(A[] array) { auto length = array.length; if (!length) return -1; auto i = uniform(0, length); auto hop = uniform(1, length); // greatest prime < 2^^31, for simplicity enum prime = 2147483477; assert (length <= prime); auto skipLength = ((prime - length) % hop) + length; auto count = 0; for (;;) { if (P(a[i])) return i; ++count; if (count == length) return -1; i += hop; if (i < length) continue; if (i < skipLength) i += hop;
(I guess this should be 'while'.)
       i -= skipLength;
    }
    return -1;
 }

 This solution is stupidly simple and I haven't tested it. But I believe
 it's truly random, since the hopping interval is arbitrary. Who knows,
 it might work.
Counterexample: [1,1,1,0,0] Your algorithm returns 0 with probability 7/20, 1 with probability 6/20 and 2 with probability 7/20. Note that there is a simple reason why your algorithm cannot work for this case: it picks one of 20 numbers at random. The resulting probability mass cannot be evenly distributed among the three elements, because 20 is not divisible by 3.
May 02
parent reply MysticZach <reachzach ggmail.com> writes:
On Tuesday, 2 May 2017 at 21:00:36 UTC, Timon Gehr wrote:
 On 02.05.2017 22:42, MysticZach wrote:
 On Tuesday, 2 May 2017 at 13:44:03 UTC, MysticZach wrote:
    for (;;) {
       if (P(a[i]))
          return i;
       ++count;
       if (count == length)
          return -1;
       i += hop;
       if (i < length)
          continue;
       if (i < skipLength)
          i += hop;
(I guess this should be 'while'.)
skipLength is determined modulo hop, thus it can't be more than one hop away.
 Counterexample: [1,1,1,0,0]

 Your algorithm returns 0 with probability 7/20, 1 with 
 probability 6/20 and 2 with probability 7/20.

 Note that there is a simple reason why your algorithm cannot 
 work for this case: it picks one of 20 numbers at random. The 
 resulting probability mass cannot be evenly distributed among 
 the three elements, because 20 is not divisible by 3.
That's true. Two points, though: If the range of error is within 1/(n*(n-1)), with array length n, it may be close enough for a given real world use case. 2. n is known to be large, which makes 1/(n*(n-1)) even smaller. You'd probably have to trade accuracy for performance. I wonder how much less performant a truly accurate algorithm would be. Also, whether there's a way to make an algorithm of this style completely accurate.
May 02
next sibling parent MysticZach <reachzach ggmail.com> writes:
On Tuesday, 2 May 2017 at 23:09:54 UTC, MysticZach wrote:
 On Tuesday, 2 May 2017 at 21:00:36 UTC, Timon Gehr wrote:
 On 02.05.2017 22:42, MysticZach wrote:
 On Tuesday, 2 May 2017 at 13:44:03 UTC, MysticZach wrote:
    for (;;) {
       if (P(a[i]))
          return i;
       ++count;
       if (count == length)
          return -1;
       i += hop;
       if (i < length)
          continue;
       if (i < skipLength)
          i += hop;
(I guess this should be 'while'.)
skipLength is determined modulo hop, thus it can't be more than one hop away.
Actually, I did botch the implementation. The hopping interval must be based on the prime rather than on n. But you could still short circuit the while loop, I think. `if (prime - n > hop) { skipLength = n + ((prime - n) % hop); }
May 03
prev sibling parent reply Timon Gehr <timon.gehr gmx.ch> writes:
On 03.05.2017 01:09, MysticZach wrote:
 Counterexample: [1,1,1,0,0]

 Your algorithm returns 0 with probability 7/20, 1 with probability
 6/20 and 2 with probability 7/20.

 Note that there is a simple reason why your algorithm cannot work for
 this case: it picks one of 20 numbers at random. The resulting
 probability mass cannot be evenly distributed among the three
 elements, because 20 is not divisible by 3.
That's true. Two points, though: If the range of error is within 1/(n*(n-1)), with array length n,
It's not though. For example, [1,1,1,0,...,0] (length 29), you get 0 and 2 each with probability 43/116, but 1 only with probability 30/116. It might be interesting to figure out how far from uniform the distribution can get.
 it may be close enough for a given
 real world use case. 2. n is known to be large, which makes 1/(n*(n-1))
 even smaller.  You'd probably have to trade accuracy for performance. I
 wonder how much less performant a truly accurate algorithm would be.
 Also, whether there's a way to make an algorithm of this style
 completely accurate.
For an accurate algorithm based on (unconditional) uniformly random sampling, the number of outcomes from sampling needs to be divisible by all numbers from 1 to n. (Because each of those could conceivably be the number of elements satisfying the predicate.)
May 04
parent reply MysticZach <reachzach ggmail.com> writes:
On Thursday, 4 May 2017 at 08:04:22 UTC, Timon Gehr wrote:
 On 03.05.2017 01:09, MysticZach wrote:
 Counterexample: [1,1,1,0,0]

 Your algorithm returns 0 with probability 7/20, 1 with 
 probability
 6/20 and 2 with probability 7/20.

 Note that there is a simple reason why your algorithm cannot 
 work for
 this case: it picks one of 20 numbers at random. The resulting
 probability mass cannot be evenly distributed among the three
 elements, because 20 is not divisible by 3.
That's true. Two points, though: If the range of error is within 1/(n*(n-1)), with array length n,
It's not though. For example, [1,1,1,0,...,0] (length 29), you get 0 and 2 each with probability 43/116, but 1 only with probability 30/116. It might be interesting to figure out how far from uniform the distribution can get.
Or how close it can get, depending on the range of intervals used. My math skill is shaky here. Maybe there's no way to deterministically jump to every element of an array with equal probability of hitting any given element satisfying a given predicate first. It sure would be cool if there were.
May 04
next sibling parent MysticZach <reachzach ggmail.com> writes:
On Thursday, 4 May 2017 at 13:19:43 UTC, MysticZach wrote:
 On Thursday, 4 May 2017 at 08:04:22 UTC, Timon Gehr wrote:
 On 03.05.2017 01:09, MysticZach wrote:
 That's true. Two points, though: If the range of error is 
 within
 1/(n*(n-1)), with array length n,
It's not though. For example, [1,1,1,0,...,0] (length 29), you get 0 and 2 each with probability 43/116, but 1 only with probability 30/116. It might be interesting to figure out how far from uniform the distribution can get.
Or how close it can get, depending on the range of intervals used. My math skill is shaky here. Maybe there's no way to deterministically jump to every element of an array with equal probability of hitting any given element satisfying a given predicate first. It sure would be cool if there were.
Within a small range of error, I mean.
May 04
prev sibling parent reply "H. S. Teoh via Digitalmars-d" <digitalmars-d puremagic.com> writes:
On Thu, May 04, 2017 at 01:19:43PM +0000, MysticZach via Digitalmars-d wrote:
[...]
 Maybe there's no way to deterministically jump to every element of an
 array with equal probability of hitting any given element satisfying a
 given predicate first. It sure would be cool if there were.
Of course there is. The random permutation approach works. In this case, the data itself is not allowed to be permuted, but equivalent semantics can be obtained by permuting an array of indices into the data. But the challenge is how efficiently you can generate a random permutation so that the cost of generating one doesn't outweigh the O(n) linear scan algorithm. I'm surprised there are no (known) incremental algorithms for generating a random permutation of 0..n that requires less than O(n) space. T -- Life is complex. It consists of real and imaginary parts. -- YHL
May 12
parent Ola Fosheim Grostad <ola.fosheim.grostad gmail.com> writes:
On Friday, 12 May 2017 at 18:43:53 UTC, H. S. Teoh wrote:
 I'm surprised there are no (known) incremental algorithms for 
 generating
 a random permutation of 0..n that requires less than O(n) space.
I've told you all you need to know...
May 12