[Boost-users] math toolkit review

From: Arnaldur Gylfason (arnaldur.gylfason_at_[hidden])
Date: 2007-04-25 11:18:48

Hi,

This is my contribution to the boost::math review.
This is a partial review of components 1) and 2) mentioned in the formal 
review announcement,
i.e. 1) the statistical distributions and 2) special functions.

?       What is your evaluation of the design?
Very good in my opinion. I think using global functions is very 
appropriate in libraries such as this one.

?       What is your evaluation of the implementation?
I did not look at the implementation but I compared it with another 
implementation 
and it compares favorably. More on this below.

?       What is your evaluation of the documentation?
Very clear. I only browsed it, but that was enough to get the feel for the 
library and start using it right away.
 
?       What is your evaluation of the potential usefulness of the 
library?
Very useful in applications that deal with anything related to statistics 
or probability.

?       Did you try to use the library?  With what compiler?  Did  you 
have any problems?
Yes I used it on a Linux box using g++ 4.0.2. No problems, except I was 
using version 1.33.0 of boost
so I had to change all use of tuple to be in namespace tuples. Took a 
short time though.

?       How much effort did you put into your evaluation? A glance? A 
quick reading? In-depth study?
I guess a quick reading plus I wrote 2 small test programs that compare 
accuracy and speed with
implementations from the R-project.

?       Are you knowledgeable about the problem domain?
I have been working in a development group in a Statistics department of a 
genetics company for several years
so I guess I'm fairly knowledgeable.

?       Do you think the library should be accepted as a Boost  library?
Yes, definitely.

The test programs are compareAccuracy.cc and compareSpeed.cc.
The programs were compiled with g++ 4.0.2 using -O3 and run on a Linux box 
with Intel Xeon 2.4 GHz.
I can give more info on this if anyone needs more info.

comparAccuracy calls several statistical distribution functions over a 
range of values,
comparing results with what similar routines from the R-project give.
The output gives the maximum abs difference and maximum of relative 
absolute difference between boost::math and R.
The results were:

norm pdf max diff:             5.55112e-17
norm pdf max rel diff:         2.20594e-16
t pdf max diff:                1.11022e-16
t pdf max rel diff:            1.17145e-15
norm cdf max diff:             1.11022e-16
norm cdf max rel diff:         2.55269e-15
t cdf max diff:                1.59317e-14
t cdf max rel diff:            3.18634e-14
F pdf max diff:                1.77636e-15
F pdf max rel diff:            7.05565e-14
F cdf max diff:                5.55112e-16
F cdf max rel diff:            1.56512e-14
chi2 pdf max diff:             5.55112e-17
chi2 pdf max rel diff:         4.26471e-14
chi2 cdf max diff:             3.46945e-16
chi2 cdf max rel diff:         6.27547e-14
exp pdf max diff:              0
exp pdf max rel diff:          0
exp cdf max diff:              0
exp cdf max rel diff:          0
gamma pdf max diff:            2.77556e-17
gamma pdf max rel diff:        3.84987e-15
gamma cdf max diff:            1.11022e-16
gamma cdf max rel diff:        1.85721e-15
lognorm pdf max diff:          1.11022e-16
lognorm pdf max rel diff:      2.21914e-16
lognorm cdf max diff:          1.11022e-16
lognorm cdf max rel diff:      1.43305e-15
cauchy pdf max diff:           0
cauchy pdf max rel diff:       0
cauchy cdf max diff:           1.11022e-16
cauchy cdf max rel diff:       2.2064e-16
weibull pdf max diff:          1.11022e-16
weibull pdf max rel diff:      2.22042e-16
weibull cdf max diff:          0
weibull cdf max rel diff:      0
norm quantile max diff:        8.88178e-16
norm quantile max rel diff:    4.00185e-16
t quantile max diff:           4.07907e-14
t quantile max rel diff:       inf
F quantile max diff:           1.11022e-15
F quantile max rel diff:       1.08231e-15
F quantile max diff:           1.11022e-15
F quantile max rel diff:       1.08231e-15
chi2 quantile max diff:        4.26326e-14
chi2 quantile max rel diff:    3.92646e-16
exp quantile max diff:         5.55112e-17
exp quantile max rel diff:     2.11093e-16
gamma quantile max diff:       3.55271e-15
gamma quantile max rel diff:   3.56265e-16
lognorm quantile max diff:     1.42109e-14
lognorm quantile max rel diff: 9.10487e-16
cauchy quantile max diff:      2.17426e-12
cauchy quantile max rel diff:  0.104537
weibull quantile max diff:     1.11022e-16
weibull quantile max rel diff: 2.0992e-16
uniform pdf max diff:          0
uniform pdf max rel diff:      0
uniform cdf max diff:          0
uniform cdf max rel diff:      0
uniform quantile max diff:     0
uniform quantile max rel diff: 0
beta pdf max diff:             1.11022e-15
beta pdf max rel diff:         3.14562e-15
beta cdf max diff:             4.44089e-16
beta cdf max rel diff:         2.9679e-15
beta quantile max diff:        2.22045e-16
beta quantile max rel diff:    8.38281e-16
binom pdf max diff:            2.08167e-17
binom pdf max rel diff:        0.00011212
binom cdf max diff:            9.4369e-16
binom cdf max rel diff:        3.46145e-14
poisson pdf max diff:          1.04083e-17
poisson pdf max rel diff:      1.05393e-13
poisson cdf max diff:          0
poisson cdf max rel diff:      0
nbinom pdf max diff:           1.38778e-17
nbinom pdf max rel diff:       1.16396e-13
nbinom cdf max diff:           8.88178e-16
nbinom cdf max rel diff:       2.37872e-15
lgamma max diff:               1.13687e-13
lgamma max rel diff:           0.211576
tgamma max diff:               7.43092e+142
tgamma max rel diff:           9.15759e-14

A few results need comment.
->t quantile max diff:           4.07907e-14
->t quantile max rel diff:       inf
Here division by 0 gives inf. max diff is quite acceptable though.
Some guards should be placed to handle 0 or near 0 but I haven't done it. 
I'm not worried about
the accuracy but anyone may feel free to remedy this.
->cauchy quantile max diff:      2.17426e-12
->cauchy quantile max rel diff:  0.104537
Here we are almost surely again dividing by something near 0. Similar 
comments as above apply.
->lgamma max diff:               1.13687e-13
->lgamma max rel diff:           0.211576
Same old story
->tgamma max diff:               7.43092e+142
->tgamma max rel diff:           9.15759e-14
At the higher end of the range tgamma gives such high numbers, that max 
diff doesn't tell us much.
max rel diff gives acceptable results though.

No quantile comparisons are given for the discrete distributions since R 
always returns integers
but boost::math gives fractional values. It depends on application what 
you want obviously
but giving results within the range (discrete in this case) is in some 
sense more natural.

All in all I am quite happy with the results.
Obviously R could be wrong but it's use is widespread in the Statistics 
community so it should be very well tested.

comparSpeed calls a few statistical distribution functions repeatedly (see 
code), timing the loop,
repeating the same thing for the R version.
The output is the time in second:

N repeats = 5000000
R norm pdf:                1.31
boost::math norm pdf:      1.1
R norm cdf:                1.49
boost::math norm cdf:      0.27
R norm quantile:           0.75
boost::math norm quantile: 1.83
R lgamma:                  4.02
boost::math lgamma:        7.4

N repeats = 1000000
R t pdf:                   0.91
boost::math t pdf:         4.85
R t cdf:                   0.17
boost::math t cdf:         0.03
R t quantile:             12.88
boost::math t quantile:   56.15
R F pdf:                   0.06
boost::math F pdf:         0.05
R F cdf:                   0.07
boost::math F cdf:         0.85
R F quantile:             18.59
boost::math F quantile:   54.54

The results vary somewhat but this would be acceptable for my 
applications.
quantile functions seem to be somewhat slower in boost::math than in R.
Maybe they're more accurate. More detailed analysis would be called for if 
people are concerned about this.

The 2 test programs are attached to this message. If they're somehow 
scrambled,
anyone can write to me and ask for them.

Last but certainly not the least, my compliments to the authors of 
boost::math ;-)

cheers

Arnaldur Gylfason,
deCODE genetics