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