Binomial Simulation

[Kenneth P. Turvey]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

I'm taking a class in Evolutionary Computation. I just turned in an
assignment that implements a simple genetic algorithm based on bit
strings. For mutation I simply check each bit in the population and
flip it with some constant probability. For a population of 10,000
individuals of 100 bits that means that I have to execute the method
Random.nextDouble() 1,000,000 times. I thought about this a bit and I
could do the same operation by getting a single value from a binomial
distribution and then selecting that number of bits at random. This
would require far fewer operations that the current brute force method.
For real problems this might not matter that much since mutation would
take a small percentage of the time spent in the algorithm, but for my
test problems mutation is the biggest contributor to run time.

So, now to the question. Is there an open source implementation of
nextBinomial() out there somewhere? I did a quick search of the web and
didn't find one. I could implement my own, but it probably wouldn't be
as fast as a widely used open source implementation and I don't want to
spend the time on something that is really a library function if I don't
have to.

Note that I'm not trying to cheat here. If I end up using it in an
assignment I will credit the author and the prof won't mind. This is a
graduate course.

Thank you.

- --
Kenneth P. Turvey <[email protected]>
http://kt.squeakydolphin.com (not much there yet)
Jabber IM: (e-mail address removed)
Phone: (314) 255-2199
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.1 (GNU/Linux)

iD8DBQFDRwbHi2ZgbrTULjoRAkxhAJsEm9g5EfvsY9P7bnE9E3czsTgsDQCg9Glz
pS2dmdxSYVmjUcfpFpUSt9A=
=vcUC
-----END PGP SIGNATURE-----

 

[zero]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

I'm taking a class in Evolutionary Computation. I just turned in an
assignment that implements a simple genetic algorithm based on bit
strings. For mutation I simply check each bit in the population and
flip it with some constant probability. For a population of 10,000
individuals of 100 bits that means that I have to execute the method
Random.nextDouble() 1,000,000 times. I thought about this a bit and I
could do the same operation by getting a single value from a binomial
distribution and then selecting that number of bits at random. This
would require far fewer operations that the current brute force method.
For real problems this might not matter that much since mutation would
take a small percentage of the time spent in the algorithm, but for my
test problems mutation is the biggest contributor to run time.

I'm not sure I understand all of this. Why would you need 1,000,000 random
values? Can't you do something like this:

for each individual
if nextDouble < probability
pick a random bit and flip it

here you would only need 10,000 random values + a fraction defined by the
probability, and not a value per bit for each of the individuals.
Typically the probability is very low, so you wouldn't have much more
random values than there are individuals in the population.

Feel free to poke fun at me if I'm missing the whole point here.

 

[Kenneth P. Turvey]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

I'm not sure I understand all of this. Why would you need 1,000,000
random values? Can't you do something like this:

for each individual
if nextDouble < probability
pick a random bit and flip it

This works as a form of mutation, but it isn't equivalent to randomly
flipping each bit with a given probability. Using this a given
individual will only have a maximum of a single bit flip.

here you would only need 10,000 random values + a fraction defined by
the
probability, and not a value per bit for each of the individuals.
Typically the probability is very low, so you wouldn't have much more
random values than there are individuals in the population.

It works, and it is certainly faster, but it isn't really the same
algorithm, and for some problems it would never help you toward a
solution. What if you really need two or three consecutive bit flips to
move you closer to the optimum. This could never happen using your
algorithm.

Feel free to poke fun at me if I'm missing the whole point here.

No, you're not. It just isn't quite what I need.

I've looked around on the web quite a bit and I haven't found anything.
I think I'm going to take the easy way out and implement it by
approximation using a normal distribution unless np and n(1-p) are less
than 5 in which case I will do it the long way, by actually pulling out
n random numbers and checking.

I hope that somebody has a better way. At least it doesn't take long to
implement.

- --
Kenneth P. Turvey <[email protected]>
http://kt.squeakydolphin.com (not much there yet)
Jabber IM: (e-mail address removed)
Phone: (314) 255-2199
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.1 (GNU/Linux)

iD8DBQFDRzLci2ZgbrTULjoRAgGjAJ40SkmhSIgyckutJt1ELZNaZn2/iwCgnp1R
Nk3Qg7z3OqK1SRmPkwn+5nw=
=S8+J
-----END PGP SIGNATURE-----

 

[Googmeister]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

I'm taking a class in Evolutionary Computation. I just turned in an
assignment that implements a simple genetic algorithm based on bit
strings. For mutation I simply check each bit in the population and
flip it with some constant probability. For a population of 10,000
individuals of 100 bits that means that I have to execute the method
Random.nextDouble() 1,000,000 times. I thought about this a bit and I
could do the same operation by getting a single value from a binomial
distribution and then selecting that number of bits at random. This
would require far fewer operations that the current brute force method.
For real problems this might not matter that much since mutation would
take a small percentage of the time spent in the algorithm, but for my
test problems mutation is the biggest contributor to run time.

Here's my understanding of your problem: given N bits, flip each
one independently with probability p, where p is some small constant.

Here's how I would approach it. Repeatedly use Geometric(p)
random variable to determine the number of bits to skip before
choosing one to flip. That is, if X_i are random geometric
random variables, then flip bits X_1, X_1 + X_2, X_1 + X_2 + X_3,
and so on. Now, each bit gets flipped with probability p,
and the memoryless property of the geometric distribution
ensures that everything is independent.

To generate a Geometric(p) variable , use
X_i = ceil(ln U / ln (1 - p)), where U is uniformly between
0 and 1. Obviously, precompute ln (1-p). Assuming p is
very small, this will be substantially more efficient than
the brute force approach.

 

[Kenneth P. Turvey]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

Here's my understanding of your problem: given N bits, flip each
one independently with probability p, where p is some small constant.

Here's how I would approach it. Repeatedly use Geometric(p)
random variable to determine the number of bits to skip before
choosing one to flip. That is, if X_i are random geometric
random variables, then flip bits X_1, X_1 + X_2, X_1 + X_2 + X_3,
and so on. Now, each bit gets flipped with probability p,
and the memoryless property of the geometric distribution
ensures that everything is independent.

That works. Now for cases where n p > 5 and n (1-p) > 5 I
can use a normal approximation and then select which bits to
flip at random. I don't know which of these two methods would
be faster. It looks like we would need the same number of
random variables roughly.

For cases where n*p < 5 or n*(1-p) < 5 then I'm looking at a very
small number of cases.

Looking at this closer, there really isn't any need to worry
about those cases I can't approximate as normal. I'm never going
to have a case where it can't be approximated. I was thinking of
n as the number of individuals in the population, but it is really
the total number of bits in the population.

Thanks for your help. I think I'll test it out the way I've
already written the nextBinomial() method and if it is too
slow, I'll try it your way.

- --
Kenneth P. Turvey <[email protected]>
http://kt.squeakydolphin.com (not much there yet)
Jabber IM: (e-mail address removed)
Phone: (314) 255-2199
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.1 (GNU/Linux)

iD8DBQFDR09Ci2ZgbrTULjoRAsi8AKDhMkTbP6wG2Ax6rBugDoRBFqP33ACgjtaY
DgfrD93vqY/bbACeGAbwnr0=
=cD39
-----END PGP SIGNATURE-----

 

[Roedy Green]

I'm taking a class in Evolutionary Computation. I just turned in an
assignment that implements a simple genetic algorithm based on bit
strings. For mutation I simply check each bit in the population and
flip it with some constant probability. For a population of 10,000
individuals of 100 bits that means that I have to execute the method
Random.nextDouble() 1,000,000 times. I thought about this a bit and I
could do the same operation by getting a single value from a binomial
distribution and then selecting that number of bits at random. This
would require far fewer operations that the current brute force method.
For real problems this might not matter that much since mutation would
take a small percentage of the time spent in the algorithm, but for my
test problems mutation is the biggest contributor to run time.

You want a bell shaped distribution, that similar to nextGaussian but
is actually a binomial. It seems to me for large n you can
approximate the binomial with a continuous function. Intuitively, I
think what you do is take the derivative of that bell shaped curve and
invert it - exchange y and x. That becomes your warping function. You
put your random 0..1 into it to get a warped random. The curve you
want is like an S lying on its side, with an arm swooping down on the
left and up on the right. It tends to move the numbers around .5 to
the mean .5 with very few near the ends mapping to numbers near 0 or
1.

If you can't solve it analytically, perhaps you can at least get the
shape of that curve through numerical methods. Then you warp your
random number to determine how many genes flipped by a table lookup
and interpolation.

 

[Roedy Green]

Thanks for your help. I think I'll test it out the way I've
already written the nextBinomial() method and if it is too
slow, I'll try it your way.

it would be interesting to hear back on a speed comparison.

 

[Chris Uppal]

For mutation I simply check each bit in the population and
flip it with some constant probability. For a population of 10,000
individuals of 100 bits that means that I have to execute the method
Random.nextDouble() 1,000,000 times. I thought about this a bit and I
could do the same operation by getting a single value from a binomial
distribution and then selecting that number of bits at random.

Do you actually need that level of statistical integrity for this application ?
It seems to me (sitting here, far away from you prof ;-) that you'd probably
get very similar results by just choosing 100*p bits (randomly) and flipping
them. Or, if that's /too/ crude, approximate it with a normal distribution ?
(Actually, it might be more interesting to try all three to see what quantitive
differences, if any, different approximations make...)

-- chris

 

[Wibble]

Kenneth P. Turvey wrote:




Do you actually need that level of statistical integrity for this application ?
It seems to me (sitting here, far away from you prof ;-) that you'd probably
get very similar results by just choosing 100*p bits (randomly) and flipping
them. Or, if that's /too/ crude, approximate it with a normal distribution ?
(Actually, it might be more interesting to try all three to see what quantitive
differences, if any, different approximations make...)

-- chris

You could treat your random as a bitmask.

hashToWidth(100, random(2^100)*probalilty)

and then just xor that with your bit population.

That would require only 10000 randoms, if you had 100 bit numbers.
Chop that into batches of 32 for 100*10000/32 randoms for 32 bit integers.

Theres some loss in statistical integrity since a given random will
hash to the same bitmap in its range, but it should be fast.

 

[Googmeister]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1



That works. Now for cases where n p > 5 and n (1-p) > 5 I
can use a normal approximation and then select which bits to
flip at random. I don't know which of these two methods would
be faster. It looks like we would need the same number of
random variables roughly.

I strongly suspect the Binomial approach (where you approximate
the Binomial with a Normal or Poisson) will be faster than the
Geometric approach I proposed since you only need a couple of
transcendental function calls instead one per bit flipped. The
advantage of the Geometric approach is that it samples from
the exact distribution (up to the randomness of random()). As
others have pointed out, such an approximation probably won't
make any difference in your application as long as np is sufficiently
large.

 

[Wibble]

You could treat your random as a bitmask.

hashToWidth(100, random(2^100)*probalilty)

and then just xor that with your bit population.

That would require only 10000 randoms, if you had 100 bit numbers.
Chop that into batches of 32 for 100*10000/32 randoms for 32 bit integers.

Theres some loss in statistical integrity since a given random will
hash to the same bitmap in its range, but it should be fast.

I couldn't resist

public final class Bits {
/* untested */
private final static int BITS_PER_WORD = 32;
private final java.util.Random r_ = new java.util.Random();
private final java.util.Random h_ = new java.util.Random();
private final float probability_;
private int bits_ = 0;
private int avail_ = 0;
public Bits(float probability) { probability_ = probability; }

public final int nextBits(int count) {
if (count<=0 || count>BITS_PER_WORD)
throw new IllegalArgumentException();
if (avail_==0) {
final int r = (int)(r_.nextInt() * probability_);
h_.setSeed(r);
bits_ = h_.nextInt();
avail_ = BITS_PER_WORD;
}
final int use = avail_<count ? avail_ : count;
final int bits = bits_ & ((1<<use)-1);
bits_ >>= use;
avail_ -= use;
return(count==use ? bits : (bits<<count)|nextBits(count-use));
}
public static void main(String argv[]) {
final int[][] population = new int[10000][];
for(int i=0;i<10000;++i) population = new int[4];
final Bits b = new Bits(0.2F/*probability*/);
for(int i=0;i<10000;++i) // process population
for(int j=0;j<4;++j) // process 100 bits
population[j] ^= b.nextBits(j<3 ? 32 : 4);
}
}

 

[Wibble]

You could treat your random as a bitmask.

hashToWidth(100, random(2^100)*probalilty)

and then just xor that with your bit population.

That would require only 10000 randoms, if you had 100 bit numbers.
Chop that into batches of 32 for 100*10000/32 randoms for 32 bit
integers.

Theres some loss in statistical integrity since a given random will
hash to the same bitmap in its range, but it should be fast.

I couldn't resist

public final class Bits {
/* untested */
private final static int BITS_PER_WORD = 32;
private final java.util.Random r_ = new java.util.Random();
private final java.util.Random h_ = new java.util.Random();
private final float probability_;
private int bits_ = 0;
private int avail_ = 0;
public Bits(float probability) { probability_ = probability; }

public final int nextBits(int count) {
if (count<=0 || count>BITS_PER_WORD)
throw new IllegalArgumentException();
if (avail_==0) {
final int r = (int)(r_.nextInt() * probability_);
h_.setSeed(r);
bits_ = h_.nextInt();
avail_ = BITS_PER_WORD;
}
final int use = avail_<count ? avail_ : count;
final int bits = bits_ & ((1<<use)-1);
bits_ >>= use;
avail_ -= use;
return(count==use ? bits : (bits<<count)|nextBits(count-use));
}
public static void main(String argv[]) {
final int[][] population = new int[10000][];
for(int i=0;i<10000;++i) population = new int[4];
final Bits b = new Bits(0.2F/*probability*/);
for(int i=0;i<10000;++i) // process population
for(int j=0;j<4;++j) // process 100 bits
population[j] ^= b.nextBits(j<3 ? 32 : 4);
}
}

oops should be
return(count==use ? bits(bits<<use/*NotCount*/)|nextBits(count-use));

 

[Kenneth P. Turvey]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

I strongly suspect the Binomial approach (where you approximate
the Binomial with a Normal or Poisson) will be faster than the
Geometric approach I proposed since you only need a couple of
transcendental function calls instead one per bit flipped. The
advantage of the Geometric approach is that it samples from
the exact distribution (up to the randomness of random()). As
others have pointed out, such an approximation probably won't
make any difference in your application as long as np is sufficiently
large.

I looked at it closer and np doesn't really have to be very big at all.
For my purposes if I use a population of over five individuals I'm
there. This will always be the case. Since others have posted code,
I'll post mine. I doubt anyone will need this, and it isn't clean, but
it passes my tests.

- --
Kenneth P. Turvey <[email protected]>
http://kt.squeakydolphin.com (not much there yet)
Jabber IM: (e-mail address removed)
Phone: (314) 255-2199

/**
* The following code has been granted to the public by the author. It
* is not covered by copyright and you have the right to do with it as
* you please.
*/

import java.util.Random;

/**
* This class implements a random number generator with some features that the
* standard Java random number generator doesn't have.
*
* @author Kenneth P. Turvey
* @version $Revision: 293 $
*/
public class ExtendedRandom extends Random {

/** Creates a new instance of ExtendedRandom */
public ExtendedRandom() {
super();
}

public ExtendedRandom(long seed) {
super(seed);
}

/**
* Returns the number of successful trials from a series of Bernouli trials
* each with the given probability.
*
* @probability the probability of each Bernouli trial.
* @trials the number of trials executed.
* @return the number of successful Bernouli trials.
*/
public int nextBinomial(int trials, double probability) {
if (trials * probability < 5.0
|| trials * (1.0 - probability) < 5.0) {
// Not enough trials for a normal approximation.

// Special case of probablity = 1/2.
// This might take less time than using nextDouble(), I
// haven't checked.
if (probability == 0.5) {
int successes = 0;
for (int count = 0; count < trials; count++) {
successes += nextBoolean() ? 1 : 0;
}
return successes;
}
else {
int successes = 0;
for (int count = 0; count < trials; count++) {
successes += super.nextDouble() < probability ? 1 : 0;
}
return successes;
}
}
else {
// We can use a normal approximation
return (int) Math.floor(
Math.sqrt(trials * probability * (1.0 - probability))
* nextGaussian()+ 0.5 + probability * trials);

}
}
}
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.1 (GNU/Linux)

iD8DBQFDR+78i2ZgbrTULjoRAt59AKCkYTmweBWjR1cbFNsoaGZ6rabEVwCfRpBY
AwSsxvcxoLTKo11YEaDOg2w=
=FhjX
-----END PGP SIGNATURE-----

 

[Chris Uppal]

// We can use a normal approximation
return (int) Math.floor(
Math.sqrt(trials * probability * (1.0 - probability))
* nextGaussian()+ 0.5 + probability * trials);

I think this can evaluate to a negative number, which is perhaps not what your
want.

-- chris

 

[Chris Uppal]

final int r = (int)(r_.nextInt() * probability_);
h_.setSeed(r);
bits_ = h_.nextInt();
avail_ = BITS_PER_WORD;

I don't fully understand what you are trying to do here but I'm fairly sure it
isn't doing what you expect.

final int r = (int)(r_.nextInt() * probability_);

nextInt() returns a value randomly chosen in the full range of signed 32-bit
integers, and then you multiply that by the probability, thus cutting down the
range. E.g if probability_ is 0.005, then the range will be roughly -10
million to +10 million.

h_.setSeed(r);

this just sets the state of the other random stream. The fact that you have
given it a restricted range of seed has no significant effect on the
subsequently generated numbers.(except to reduce the number of possible
different /streams/ of numbers that _h can produce).

bits_ = h_.nextInt();

this will, again, return a value randomly chosen in the full range of signed
32-bit integers. So each bit has the same probability (1/2) of being set. In
the case where probability_ is 0.005, there are only 20 million different seed
values that you can have used, and so there are only 20 million different
possible sequences of 32 bits, and hence only that many possible values for
bits_ (instead of roughly 4 billion -- which isn't ideal, but is probably not
terribly important here).

I /think/ that what you are intending to do is generate an integer such that
each bit has a 'probability_' chance of being set, but I don't know how to do
that except by generating each bit separately, or using one of the
approximations we've been discussing.

-- chris

 

[Kenneth P. Turvey]

-----BEGIN PGP SIGNED MESSAGE-----
Hash: SHA1

I think this can evaluate to a negative number, which is perhaps not
what your want.

You're right. It is unlikely, but possible. I'll have to look at this
again.

Thanks for the tip.

I added the check to my unit tests and it immediately failed.
Apparently it isn't as unlikely as I would have thought intuitively.

The other problem is that the result can be larger than the number of
trials, which we really don't want either. I'll have to play with the
math a little bit before I decide how to handle this. I'm not sure
simple truncation of the tails (returning either 0 or the number of
trials) is really a good way to handle this.

You know, I haven't really played with any statistics in years and I'm
enjoying it. I'm not sure what that says about me.

- --
Kenneth P. Turvey <[email protected]>
http://kt.squeakydolphin.com (not much there yet)
Jabber IM: (e-mail address removed)
Phone: (314) 255-2199
-----BEGIN PGP SIGNATURE-----
Version: GnuPG v1.4.1 (GNU/Linux)

iD8DBQFDSYp3i2ZgbrTULjoRAve/AKCGXP00KZo+daz4C94znTK8DQeEeQCfYTzn
fukhD8FV0kh+fx6l1UxRtTg=
=iwij
-----END PGP SIGNATURE-----

 

[Chris Uppal]

You know, I haven't really played with any statistics in years and I'm
enjoying it. I'm not sure what that says about me.

<chuckle/>

I always had great difficulty with the subject myself (so beware of my
reasoning in what follows !).

The other problem is that the result can be larger than the number of
trials, which we really don't want either. I'll have to play with the
math a little bit before I decide how to handle this. I'm not sure
simple truncation of the tails (returning either 0 or the number of
trials) is really a good way to handle this.

A way of thinking about this that appeals to me is to anthropomorphise the
situation a bit (as, of course, we OO programmers always do). Imagine that one
of your creatures (lineages, bitmaps, whatever you want to call 'em) is --
unlike me -- a minimally competent statistician. It knows that God hands out
mutations at random, but that's /all/ it knows, except how many times it has
been zapped in each generation. So if in one generation it finds itself with
three mutations, it knows that God must have decreed that there be at least
three mutations in the population as a whole, but it doesn't now any more than
that. In particular it doesn't know how God choose the total number (call it
M), nor how M is distributed. M could, for instance, be a constant if God
didn't care for mucking around with complex code.

The question is, how long would it take the statistician to be able to tell
that God wasn't choosing M according to a pure binomial distribution ? And
would it be able to amass that much data in one run of your simulation ?
Because if it can't tell, then the difference can't be effecting its behaviour
(to a statistically significant degree) -- and its behaviour is all that
matters.

If God were using a constant value for M, then the tail of the distribution it
saw (of the number of hits it took per generation) would drop to zero
anomalously quickly (it would, for instance, /never/ take M+1 hits). My guess
is that the effect would be pretty subtle, and hard to detect in the limited
amount of data it could collect. I'm sure you will enjoy working out the exact
details ;-) If God were using a normal distribution for M (with more-or-less
arbitrary truncation), then presumably the divergence would be even harder for
the creature to measure, and would require more data (and a longer run) to
detect. If God went even further (as might not be a bad idea) and simply
pre-populated a largish array with integer values drawn from a binomial
distribution, and then chose an arbitrary[*] value from that array to use as M
for each generation, then it would -- I think -- be still harder to detect.

([*] cycle through 'em, choose one at random, I doubt whether it would make
much difference...)

-- chris