|
Home | Libraries | People | FAQ | More |
Discrete distributions present us with a problem when calculating the quantile: we are starting from a continuous real-valued variable - the probability - but the result (the value of the random variable) should really be discrete.
Consider for example a Binomial distribution, with a sample size of 50, and a success fraction of 0.5. There are a variety of ways we can plot a discrete distribution, but if we plot the PDF as a step-function then it looks something like this:
Now lets suppose that the user asks for a the quantile that corresponds to a probability of 0.05, if we zoom in on the CDF for that region here's what we see:
As can be seen there is no random variable that corresponds to a probability of exactly 0.05, so we're left with two choices as shown in the figure:
In fact there's actually a third choice as well: we could "pretend" that the distribution was continuous and return a real valued result: in this case we would calculate a result of approximately 18.701 (this accurately reflects the fact that the result is nearer to 19 than 18).
By using policies we can offer any of the above as options, but that still leaves the question: What is actually the right thing to do?
And in particular: What policy should we use by default?
Calculating an integer result is often much faster than calculating a real-valued
result (in our tests it was up to 20 times faster), and the random variable
can only take integer values anyway, so the default is an integer. Which
integer? The usual definition of the quantile of a discrete distribution
is the smallest value k for which cdf(dist,
k) >= p.
It is what R, SciPy, Julia and Mathematica return, it is the empirical quantile
of a large sample, and it is the default here: the integer_round_up
policy. In the example above the result is 19, since:
cdf(binomial(50, 0.5), 18) < 0.05 <= cdf(binomial(50, 0.5), 19)
Upper quantiles are best found from the complement: quantile(complement(dist,
q))
is the smallest k for which cdf(complement(dist,
k)) <= q,
that is, for which the probability of a value greater than k
is at most q. This is quantile(dist,
1 - q) without
the rounding error in 1 -
q, which matters when q
is small. For example:
quantile(complement(binomial(50, 0.5), 0.05));
returns 31, since:
cdf(complement(binomial(50, 0.5), 30)) > 0.05 >= cdf(complement(binomial(50, 0.5), 31))
Together these give a two-sided interval. If L and U are the lower and upper quantiles, then less than 5% of the probability lies below L and no more than 5% above U, so at least 90% lies in the central region:
L <= X && X <= U
and we would reject the null hypothesis if:
X < L || X > U
The moral here is to always be very careful with your comparisons when dealing with a discrete distribution, and if in doubt, base your comparisons on CDF's instead.
The integer_round_outwards
policy, which you can select by defining BOOST_MATH_DISCRETE_QUANTILE_POLICY
as integer_round_outwards
(see Using Macros to
Change the Policy Defaults) or by using the policy directly, rounds
lower quantiles (probability less than 0.5) down and upper quantiles (probability
greater than 0.5) up. The reasoning is that if someone asks for a quantile
at 0.05, then we can normally assume that they are asking for at least 95% of the probability to the right of the
value chosen, and no more than 5% of the
probability to the left of the value chosen, counting the value
chosen as part of the left tail.
In the above binomial example the 0.05 quantile is then rounded down to 18. The converse applies to upper-quantiles: they are rounded up, so that at least the requested probability is to the left of the value returned, and no more than 1 - the requested probability is to the right of the value returned.
For example, taking our 50 sample binomial distribution with a success fraction of 0.5, if we wanted a two sided 90% confidence interval, then we would ask for the 0.05 and 0.95 quantiles with the results rounded outwards so that at least 90% of the probability is in the central area:
With this policy:
quantile(binomial(50, 0.5), 0.05);
returns 18, and there is no x greater than 18 for which:
cdf(binomial(50, 0.5), x) <= 0.05;
However:
quantile(binomial(50, 0.5), 0.95);
returns 31, and there is no x less than 31 for which:
cdf(binomial(50, 0.5), x) >= 0.95;
We might naively expect that for this symmetrical distribution the result would be 32 (since 32 = 50 - 18), but we need to remember that the cdf of the binomial is inclusive of the random variable. So while the left tail area includes the quantile returned, the right tail area always excludes an upper quantile value: since that "belongs" to the central area. Therefore if U and L are the upper and lower quantiles respectively, then a random variable X is in the tail area - where we would reject the null hypothesis if:
X <= L || X > U
And the a variable X is inside the central region if:
L < X <= U
Unless the lower probability is exactly the cdf at some integer, this is the same region as the default gives: only the lower end is reported one lower. Note also that the result jumps at a probability of 0.5, where the rounding changes direction, so (for example) the quantiles of uniformly spaced probabilities are not spread evenly around the median.
As you would expect from a section on policies, you won't be surprised to know that other rounding options are available:
As described above: lower quantiles are rounded down (probability < 0.5), and upper quantiles (probability > 0.5) are rounded up.
This gives no more than the requested probability in the tails, and at least the requested probability in the central area.
This is the exact opposite of integer_round_outwards:
lower quantiles are rounded up (probability < 0.5), and upper quantiles
(probability > 0.5) are rounded down.
This gives at least the requested probability in the tails, and no more than the requested probability in the central area.
This policy will always round the result down no matter whether it
is an upper or lower quantile, giving the largest k
with cdf(dist, k) <= p
(or cdf(complement(dist, k)) >= q)
This is the default policy: it will always round the result up no matter
whether it is an upper or lower quantile, giving the smallest k
with cdf(dist, k) >= p
(or cdf(complement(dist, k)) <= q)
This policy will always round the result to the nearest integer no matter whether it is an upper or lower quantile
This policy will return a real valued result for the quantile of a discrete distribution: this is generally much slower than finding an integer result but does allow for more sophisticated rounding policies.
To understand how the rounding policies for the discrete distributions can be used, we'll use the 50-sample binomial distribution with a success fraction of 0.5 once again, and calculate all the possible quantiles at 0.05 and 0.95.
Begin by including the needed headers (and some using statements for conciseness):
#include <iostream> using std::cout; using std::endl; using std::left; using std::fixed; using std::right; using std::scientific; #include <iomanip> using std::setw; using std::setprecision; #include <boost/math/distributions/binomial.hpp>
Next we'll bring the needed declarations into scope, and define distribution types for all the available rounding policies:
// Avoid // using namespace std; // and // using namespace boost::math; // to avoid potential ambiguity of names, like binomial. // using namespace boost::math::policies; is small risk, but // the necessary items are brought into scope thus: using boost::math::binomial_distribution; using boost::math::policies::policy; using boost::math::policies::discrete_quantile; using boost::math::policies::integer_round_outwards; using boost::math::policies::integer_round_down; using boost::math::policies::integer_round_up; using boost::math::policies::integer_round_nearest; using boost::math::policies::integer_round_inwards; using boost::math::policies::real; using boost::math::binomial_distribution; // Not std::binomial_distribution. typedef binomial_distribution< double, policy<discrete_quantile<integer_round_outwards> > > binom_round_outwards; typedef binomial_distribution< double, policy<discrete_quantile<integer_round_inwards> > > binom_round_inwards; typedef binomial_distribution< double, policy<discrete_quantile<integer_round_down> > > binom_round_down; typedef binomial_distribution< double, policy<discrete_quantile<integer_round_up> > > binom_round_up; typedef binomial_distribution< double, policy<discrete_quantile<integer_round_nearest> > > binom_round_nearest; typedef binomial_distribution< double, policy<discrete_quantile<real> > > binom_real_quantile;
Now let's set to work calling those quantiles:
int main() { cout << "Testing rounding policies for a 50 sample binomial distribution,\n" "with a success fraction of 0.5.\n\n" "Lower quantiles are calculated at p = 0.05\n\n" "Upper quantiles at p = 0.95.\n\n"; cout << setw(25) << right << "Policy"<< setw(18) << right << "Lower Quantile" << setw(18) << right << "Upper Quantile" << endl; // Test integer_round_outwards: cout << setw(25) << right << "integer_round_outwards" << setw(18) << right << quantile(binom_round_outwards(50, 0.5), 0.05) << setw(18) << right << quantile(binom_round_outwards(50, 0.5), 0.95) << endl; // Test integer_round_inwards: cout << setw(25) << right << "integer_round_inwards" << setw(18) << right << quantile(binom_round_inwards(50, 0.5), 0.05) << setw(18) << right << quantile(binom_round_inwards(50, 0.5), 0.95) << endl; // Test integer_round_down: cout << setw(25) << right << "integer_round_down" << setw(18) << right << quantile(binom_round_down(50, 0.5), 0.05) << setw(18) << right << quantile(binom_round_down(50, 0.5), 0.95) << endl; // Test integer_round_up: cout << setw(25) << right << "integer_round_up" << setw(18) << right << quantile(binom_round_up(50, 0.5), 0.05) << setw(18) << right << quantile(binom_round_up(50, 0.5), 0.95) << endl; // Test integer_round_nearest: cout << setw(25) << right << "integer_round_nearest" << setw(18) << right << quantile(binom_round_nearest(50, 0.5), 0.05) << setw(18) << right << quantile(binom_round_nearest(50, 0.5), 0.95) << endl; // Test real: cout << setw(25) << right << "real" << setw(18) << right << quantile(binom_real_quantile(50, 0.5), 0.05) << setw(18) << right << quantile(binom_real_quantile(50, 0.5), 0.95) << endl; } // int main()
Which produces the program output:
policy_eg_10.vcxproj -> J:\Cpp\MathToolkit\test\Math_test\Release\policy_eg_10.exe
Testing rounding policies for a 50 sample binomial distribution,
with a success fraction of 0.5.
Lower quantiles are calculated at p = 0.05
Upper quantiles at p = 0.95.
Policy Lower Quantile Upper Quantile
integer_round_outwards 18 31
integer_round_inwards 19 30
integer_round_down 18 30
integer_round_up 19 31
integer_round_nearest 19 30
real 18.701 30.299