[random] A much faster Poisson variate generation

Hello, I have modified the poisson_distribution class to use an updated version of Ahrens and Dieter's PD algorithm for generating Poisson variates, as described in "Computer Generation of Poisson Deviates from Modified Normal Distributions." Is this allowed in Boost, and is there interest? The current algorithm is O(mean), while PD is bounded time. In addition, the current implementation fails for mean >> 750. PD does, however, use a lookup table and thus has a non-trivial constructor and needs 36 * sizeof(RealType) + 1 * sizeof(IntType) bytes per poisson_distribution object. Please let me know if this is desired. Ross

-----Original Message----- From: boost-bounces@lists.boost.org [mailto:boost-bounces@lists.boost.org] On Behalf Of Ross Levine Sent: 30 December 2008 19:42 To: boost@lists.boost.org Subject: [boost] [random] A much faster Poisson variate generation
Hello, I have modified the poisson_distribution class to use an updated version of Ahrens and Dieter's PD algorithm for generating Poisson variates, as described in "Computer Generation of Poisson Deviates from Modified Normal Distributions." Is this allowed in Boost, and is there interest? The current algorithm is O(mean), while PD is bounded time. In addition, the current implementation fails for mean >> 750. PD does, however, use a lookup table and thus has a non-trivial constructor and needs 36 * sizeof(RealType) + 1 * sizeof(IntType) bytes per poisson_distribution object. Please let me know if this is desired.
Is this http://doi.acm.org/10.1145/355993.355997 aka Computer Generation of Poisson Deviates from Modified Normal Distributions Full text PdfPdf (791 KB) Source ACM Transactions on Mathematical Software (TOMS) archive Volume 8 , Issue 2 (June 1982) table of contents Pages: 163 - 179 Year of Publication: 1982 ISSN:0098-3500 Authors J. H. Ahrens Department of Mathematics, University of Kiel, D 2300 Kiel, Olshausenstrasse 40-60, West Germany U. Dieter Institute of Statistics, Technical University, A 8010 Graz, Hamerlinggasse 6, Austria Publisher ACM New York, NY, USA Bibliometrics Downloads (6 Weeks): 34, Downloads (12 Months): 280, Citation Count: 5 (This requires a subscription to ACM) Or is it a later improved version? Sounds useful, but I expect can you provide some tests to show equivalence and enhancedness? Paul --- Paul A. Bristow Prizet Farmhouse Kendal, UK LA8 8AB +44 1539 561830, mobile +44 7714330204 pbristow@hetp.u-net.com

That is the article in question. The original algorithm was in a pseudo-assembly language. I have updated it to use binary search and a few other optimizations. I have tested the class extensively. It is about 15 times faster than boost::poisson_distribution for _mean=500. At small (<10) means both are too fast to measure. As expected, the mean and standard distribution of samples from this poisson distribution are equal to _mean. I'm not sure exactly how to post code here. I'm afraid the spaces will be deleted or something. I'll send it in my next message. On Wed, Dec 31, 2008 at 1:36 PM, Paul A. Bristow <pbristow@hetp.u-net.com> wrote:
-----Original Message----- From: boost-bounces@lists.boost.org [mailto:boost-bounces@lists.boost.org] On Behalf Of Ross Levine Sent: 30 December 2008 19:42 To: boost@lists.boost.org Subject: [boost] [random] A much faster Poisson variate generation
Hello, I have modified the poisson_distribution class to use an updated version of Ahrens and Dieter's PD algorithm for generating Poisson variates, as described in "Computer Generation of Poisson Deviates from Modified Normal Distributions." Is this allowed in Boost, and is there interest? The current algorithm is O(mean), while PD is bounded time. In addition, the current implementation fails for mean >> 750. PD does, however, use a lookup table and thus has a non-trivial constructor and needs 36 * sizeof(RealType) + 1 * sizeof(IntType) bytes per poisson_distribution object. Please let me know if this is desired.
Is this http://doi.acm.org/10.1145/355993.355997
aka
Computer Generation of Poisson Deviates from Modified Normal Distributions Full text PdfPdf (791 KB) Source ACM Transactions on Mathematical Software (TOMS) archive Volume 8 , Issue 2 (June 1982) table of contents Pages: 163 - 179 Year of Publication: 1982 ISSN:0098-3500 Authors J. H. Ahrens Department of Mathematics, University of Kiel, D 2300 Kiel, Olshausenstrasse 40-60, West Germany U. Dieter Institute of Statistics, Technical University, A 8010 Graz, Hamerlinggasse 6, Austria Publisher ACM New York, NY, USA Bibliometrics Downloads (6 Weeks): 34, Downloads (12 Months): 280, Citation Count: 5
(This requires a subscription to ACM)
Or is it a later improved version?
Sounds useful, but I expect can you provide some tests to show equivalence and enhancedness?
Paul
--- Paul A. Bristow Prizet Farmhouse Kendal, UK LA8 8AB +44 1539 561830, mobile +44 7714330204 pbristow@hetp.u-net.com
_______________________________________________ Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost

Silly me, I could always just link to it... http://dl.getdropbox.com/u/416909/poisson_distribution.hpp And the source file I use to test... http://dl.getdropbox.com/u/416909/test_pd.cpp Of special note is that the boost version of the poisson distribution does not work above mean=745 on my machine, it seems to encounter a rounding error. On Wed, Dec 31, 2008 at 3:28 PM, Ross Levine <ross.levine@uky.edu> wrote:
That is the article in question. The original algorithm was in a pseudo-assembly language. I have updated it to use binary search and a few other optimizations. I have tested the class extensively. It is about 15 times faster than boost::poisson_distribution for _mean=500. At small (<10) means both are too fast to measure.
As expected, the mean and standard distribution of samples from this poisson distribution are equal to _mean. I'm not sure exactly how to post code here. I'm afraid the spaces will be deleted or something. I'll send it in my next message.
On Wed, Dec 31, 2008 at 1:36 PM, Paul A. Bristow <pbristow@hetp.u-net.com> wrote:
-----Original Message----- From: boost-bounces@lists.boost.org [mailto:boost-bounces@lists.boost.org] On Behalf Of Ross Levine Sent: 30 December 2008 19:42 To: boost@lists.boost.org Subject: [boost] [random] A much faster Poisson variate generation
Hello, I have modified the poisson_distribution class to use an updated version of Ahrens and Dieter's PD algorithm for generating Poisson variates, as described in "Computer Generation of Poisson Deviates from Modified Normal Distributions." Is this allowed in Boost, and is there interest? The current algorithm is O(mean), while PD is bounded time. In addition, the current implementation fails for mean >> 750. PD does, however, use a lookup table and thus has a non-trivial constructor and needs 36 * sizeof(RealType) + 1 * sizeof(IntType) bytes per poisson_distribution object. Please let me know if this is desired.
Is this http://doi.acm.org/10.1145/355993.355997
aka
Computer Generation of Poisson Deviates from Modified Normal Distributions Full text PdfPdf (791 KB) Source ACM Transactions on Mathematical Software (TOMS) archive Volume 8 , Issue 2 (June 1982) table of contents Pages: 163 - 179 Year of Publication: 1982 ISSN:0098-3500 Authors J. H. Ahrens Department of Mathematics, University of Kiel, D 2300 Kiel, Olshausenstrasse 40-60, West Germany U. Dieter Institute of Statistics, Technical University, A 8010 Graz, Hamerlinggasse 6, Austria Publisher ACM New York, NY, USA Bibliometrics Downloads (6 Weeks): 34, Downloads (12 Months): 280, Citation Count: 5
(This requires a subscription to ACM)
Or is it a later improved version?
Sounds useful, but I expect can you provide some tests to show equivalence and enhancedness?
Paul
--- Paul A. Bristow Prizet Farmhouse Kendal, UK LA8 8AB +44 1539 561830, mobile +44 7714330204 pbristow@hetp.u-net.com
_______________________________________________ Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost

Your test promises equivalence (and much improvement) over the existing code, but are you able to run the existing test? \boost_1_37_0\libs\random\random_test.cpp Which tests with instantiate_dist(urng, "poisson_distribution", boost::poisson_distribution<>(1)); Are there any limits that can be relaxed now that you have 'debugged' the 'biggish mean' case? Are your changes likely to cause any problems to existing users? (breaking existing systems is unpopular ;-) If this looks OK, I suspect that someone can just copy this into trunk and check if all is well? John Maddock? Paul PS The existing authors copyrights in the header should be retained, and yours added of course, with some notes on the changes made (since they are quite major in implementation). And some notes on your changes need to be added to the docs at /libs/random/index.html Would you like to suggest working to go in the History and Acknowledgements section. I am not clear what tool was used to generate the html, but I suspect we could cheat and edited by hand for such a small change.
-----Original Message----- From: boost-bounces@lists.boost.org [mailto:boost-bounces@lists.boost.org] On Behalf Of Ross Levine Sent: 31 December 2008 20:43 To: boost@lists.boost.org Subject: Re: [boost] [random] A much faster Poisson variate generation
Silly me, I could always just link to it...
http://dl.getdropbox.com/u/416909/poisson_distribution.hpp
And the source file I use to test...
http://dl.getdropbox.com/u/416909/test_pd.cpp
Of special note is that the boost version of the poisson distribution does not work above mean=745 on my machine, it seems to encounter a rounding error.
On Wed, Dec 31, 2008 at 3:28 PM, Ross Levine <ross.levine@uky.edu> wrote:
That is the article in question. The original algorithm was in a pseudo-assembly language. I have updated it to use binary search and a few other optimizations. I have tested the class extensively. It is about 15 times faster than boost::poisson_distribution for _mean=500. At small (<10) means both are too fast to measure.
As expected, the mean and standard distribution of samples from this poisson distribution are equal to _mean. I'm not sure exactly how to post code here. I'm afraid the spaces will be deleted or something. I'll send it in my next message.
On Wed, Dec 31, 2008 at 1:36 PM, Paul A. Bristow <pbristow@hetp.u-net.com> wrote:
-----Original Message----- From: boost-bounces@lists.boost.org
[mailto:boost-bounces@lists.boost.org]
On
Behalf Of Ross Levine Sent: 30 December 2008 19:42 To: boost@lists.boost.org Subject: [boost] [random] A much faster Poisson variate generation
Hello, I have modified the poisson_distribution class to use an updated version of Ahrens and Dieter's PD algorithm for generating Poisson variates, as described in "Computer Generation of Poisson Deviates from Modified Normal Distributions." Is this allowed in Boost, and is there interest? The current algorithm is O(mean), while PD is bounded time. In addition, the current implementation fails for mean >> 750. PD does, however, use a lookup table and thus has a non-trivial constructor and needs 36 * sizeof(RealType) + 1 * sizeof(IntType) bytes per poisson_distribution object. Please let me know if this is desired.
Is this http://doi.acm.org/10.1145/355993.355997
aka
Computer Generation of Poisson Deviates from Modified Normal Distributions Full text PdfPdf (791 KB) Source ACM Transactions on Mathematical Software (TOMS) archive Volume 8 , Issue 2 (June 1982) table of contents Pages: 163 - 179 Year of Publication: 1982 ISSN:0098-3500 Authors J. H. Ahrens Department of Mathematics, University of Kiel, D 2300 Kiel, Olshausenstrasse 40-60, West Germany U. Dieter Institute of Statistics, Technical University, A 8010 Graz, Hamerlinggasse 6, Austria Publisher ACM New York, NY, USA Bibliometrics Downloads (6 Weeks): 34, Downloads (12 Months): 280, Citation Count: 5
(This requires a subscription to ACM)
Or is it a later improved version?
Sounds useful, but I expect can you provide some tests to show equivalence and enhancedness?
Paul
--- Paul A. Bristow Prizet Farmhouse Kendal, UK LA8 8AB +44 1539 561830, mobile +44 7714330204 pbristow@hetp.u-net.com
_______________________________________________ Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost
_______________________________________________ Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost

-----Original Message----- From: boost-bounces@lists.boost.org [mailto:boost-bounces@lists.boost.org] On Behalf Of Paul A. Bristow Sent: 05 January 2009 15:25 To: boost@lists.boost.org Subject: Re: [boost] [random] A much faster Poisson variate generation
I have also just realized that your new version using MPL may have implications for 'obselete' compilers. It would be best if the old simpler version still worked for these compilers. Perhaps you have already given this some thought - I see some tests for MSVC <=1300 ? At least, it will require altering the Boost 'expected to fail' list. Paul --- Paul A. Bristow Prizet Farmhouse Kendal, UK LA8 8AB +44 1539 561830, mobile +44 7714330204 pbristow@hetp.u-net.com

I have removed the MPL: the code was there to use an int instead of a long for the factorial lookup table if the size of an int was large enough. Ultimately, the rare case in which this would be useful (when 8 <= sizeof(int) < sizeof(long)) is not worth code that won't compile on some compilers. The class now has an empty reset() member function for compatibility, and passes random_test.cpp's test_main function as well, and can serve as a drop-in replacement of the existing poisson_distribution.hpp file. Notably, operator<< and operator>> are compatible with the old version, streaming only the mean. I am unaware of how to officially submit this header, and where to change documentation. Thanks, Ross On Tue, Jan 6, 2009 at 5:19 AM, Paul A. Bristow <pbristow@hetp.u-net.com> wrote:
-----Original Message----- From: boost-bounces@lists.boost.org [mailto:boost-bounces@lists.boost.org] On Behalf Of Paul A. Bristow Sent: 05 January 2009 15:25 To: boost@lists.boost.org Subject: Re: [boost] [random] A much faster Poisson variate generation
I have also just realized that your new version using MPL may have implications for 'obselete' compilers.
It would be best if the old simpler version still worked for these compilers.
Perhaps you have already given this some thought - I see some tests for MSVC <=1300 ?
At least, it will require altering the Boost 'expected to fail' list.
Paul
--- Paul A. Bristow Prizet Farmhouse Kendal, UK LA8 8AB +44 1539 561830, mobile +44 7714330204 pbristow@hetp.u-net.com
_______________________________________________ Unsubscribe & other changes: http://lists.boost.org/mailman/listinfo.cgi/boost

-----Original Message----- From: boost-bounces@lists.boost.org [mailto:boost-bounces@lists.boost.org] On Behalf Of Ross Levine Sent: 14 January 2009 20:13 To: boost@lists.boost.org Subject: Re: [boost] [random] A much faster Poisson variate generation
I have removed the MPL: the code was there to use an int instead of a long for the factorial lookup table if the size of an int was large enough. Ultimately, the rare case in which this would be useful (when 8 <= sizeof(int) < sizeof(long)) is not worth code that won't compile on some compilers. The class now has an empty reset() member function for compatibility, and passes random_test.cpp's test_main function as well, and can serve as a drop-in replacement of the existing poisson_distribution.hpp file. Notably, operator<< and operator>> are compatible with the old version, streaming only the mean.
Sounds good to me, but
I am unaware of how to officially submit this header, and where to change documentation.
I think you should contact the original author if possible - direct as your note hasn't produced a response. Ask again on the list if you fail to contact him. Paul --- Paul A. Bristow Prizet Farmhouse Kendal, UK LA8 8AB +44 1539 561830, mobile +44 7714330204 pbristow@hetp.u-net.com
participants (2)
-
Paul A. Bristow
-
Ross Levine