recombination chunk of length infinity

[GoogleCodeExporter]

In macs, it is possible to draw a distance of +infinity to the next 
recombination breakpoint, and hence for there to be no recombination for most 
or all of an entire sequence (which happened to me in one simulation).

My apologies, but this problem is difficult to reproduce, since I have slightly 
modified my version of macs. However, this error occurs with probability 2^-32 
per recombination chunk (see below), so should also affect other users. I can 
send you my modified macs if you wish -- there are only ~5 lines of code 
modified, but they do affect the order of calls to the random number generator.

The cause of the error is as follows. The distance to the next recombination 
breakpoint is gotten by drawing from an exponential distribution, which in turn 
is gotten by drawing u ~ Uniform[0,1), and transforming to log(u) / Lambda. 
However, the random number generator that macs uses has a 2^-32 probability of 
returning exactly 0, which will cause log(u) / Lambda = Infinity.

Replacing log(u) / Lambda with log(1 - u) / Lambda will fix this problem, since 
the random number generator boost:uniform_01 will never return exactly 1. This 
change can be made at the end of simulator.h.

Original issue reported on code.google.com by jackk...@gmail.com on 8 May 2013 at 5:18

[GoogleCodeExporter]

-Jack Kamm
Ph.d. student
Song group
UC Berkeley

Original comment by jackk...@gmail.com on 8 May 2013 at 5:19

[GoogleCodeExporter]

Jack, 

Great idea.  Thanks for spotting this.  I will make this change to the header 
as you suggested.

Gary

Original comment by gche...@gmail.com on 8 May 2013 at 5:20

• Changed state: Started

[GoogleCodeExporter]

Jack, let me know if this version works for you.  It's now in SVN.

Original comment by gche...@gmail.com on 8 May 2013 at 5:28

[GoogleCodeExporter]

Hi Gary,

Yes, this works. Thanks for addressing this! I was quite surprised to encounter 
such a low probability event, but given the number of simulations I was doing 
it wound up being in the realm of possibility.

--Jack

Original comment by jackk...@gmail.com on 8 May 2013 at 7:49

[GoogleCodeExporter]

Original comment by gche...@gmail.com on 7 Jun 2013 at 4:24

• Changed state: Fixed