...one of the most highly
regarded and expertly designed C++ library projects in the
world.
— Herb Sutter and Andrei
Alexandrescu, C++
Coding Standards
First we need some includes to access the Normal Distribution, the algorithms to find scale (and some std output of course).
#include <boost/math/distributions/normal.hpp> // for normal_distribution using boost::math::normal; // typedef provides default type is double. #include <boost/math/distributions/find_scale.hpp> using boost::math::find_scale; using boost::math::complement; // Needed if you want to use the complement version. using boost::math::policies::policy; // Needed to specify the error handling policy. #include <iostream> using std::cout; using std::endl; #include <iomanip> using std::setw; using std::setprecision; #include <limits> using std::numeric_limits;
For this example, we will use the standard Normal Distribution, with location (mean) zero and standard deviation (scale) unity. Conveniently, this is also the default for this implementation's constructor.
normal N01; // Default 'standard' normal distribution with zero mean double sd = 1.; // and standard deviation is 1.
Suppose we want to find a different normal distribution with standard deviation so that only fraction p (here 0.001 or 0.1%) are below a certain chosen limit (here -2. standard deviations).
double z = -2.; // z to give prob p double p = 0.001; // only 0.1% below z = -2 cout << "Normal distribution with mean = " << N01.location() // aka N01.mean() << ", standard deviation " << N01.scale() // aka N01.standard_deviation() << ", has " << "fraction <= " << z << ", p = " << cdf(N01, z) << endl; cout << "Normal distribution with mean = " << N01.location() << ", standard deviation " << N01.scale() << ", has " << "fraction > " << z << ", p = " << cdf(complement(N01, z)) << endl; // Note: uses complement.
Normal distribution with mean = 0 has fraction <= -2, p = 0.0227501 Normal distribution with mean = 0 has fraction > -2, p = 0.97725
Noting that p = 0.02 instead of our target of 0.001, we can now use
find_scale
to give
a new standard deviation.
double l = N01.location(); double s = find_scale<normal>(z, p, l); cout << "scale (standard deviation) = " << s << endl;
that outputs:
scale (standard deviation) = 0.647201
showing that we need to reduce the standard deviation from 1. to 0.65.
Then we can check that we have achieved our objective by constructing a new distribution with the new standard deviation (but same zero mean):
normal np001pc(N01.location(), s);
And re-calculating the fraction below (and above) our chosen limit.
cout << "Normal distribution with mean = " << l << " has " << "fraction <= " << z << ", p = " << cdf(np001pc, z) << endl; cout << "Normal distribution with mean = " << l << " has " << "fraction > " << z << ", p = " << cdf(complement(np001pc, z)) << endl;
Normal distribution with mean = 0 has fraction <= -2, p = 0.001 Normal distribution with mean = 0 has fraction > -2, p = 0.999
We can also control the policy for handling various errors. For example, we can define a new (possibly unwise) policy to ignore domain errors ('bad' arguments).
Unless we are using the boost::math namespace, we will need:
using boost::math::policies::policy; using boost::math::policies::domain_error; using boost::math::policies::ignore_error;
Using a typedef is convenient, especially if it is re-used, although it is not required, as the various examples below show.
typedef policy<domain_error<ignore_error> > ignore_domain_policy; // find_scale with new policy, using typedef. l = find_scale<normal>(z, p, l, ignore_domain_policy()); // Default policy policy<>, needs using boost::math::policies::policy; l = find_scale<normal>(z, p, l, policy<>()); // Default policy, fully specified. l = find_scale<normal>(z, p, l, boost::math::policies::policy<>()); // New policy, without typedef. l = find_scale<normal>(z, p, l, policy<domain_error<ignore_error> >());
If we want to express a probability, say 0.999, that is a complement,
1 -
p
we should not even think
of writing find_scale<normal>(z, 1 - p, l)
, but instead,
use the complements version.
z = -2.; double q = 0.999; // = 1 - p; // complement of 0.001. sd = find_scale<normal>(complement(z, q, l)); normal np95pc(l, sd); // Same standard_deviation (scale) but with mean(scale) shifted cout << "Normal distribution with mean = " << l << " has " << "fraction <= " << z << " = " << cdf(np95pc, z) << endl; cout << "Normal distribution with mean = " << l << " has " << "fraction > " << z << " = " << cdf(complement(np95pc, z)) << endl;
Sadly, it is all too easy to get probabilities the wrong way round, when you may get a warning like this:
Message from thrown exception was: Error in function boost::math::find_scale<Dist, Policy>(complement(double, double, double, Policy)): Computed scale (-0.48043523852179076) is <= 0! Was the complement intended?
The default error handling policy is to throw an exception with this message, but if you chose a policy to ignore the error, the (impossible) negative scale is quietly returned.
See find_scale_example.cpp for full source code: the program output looks like this:
Example: Find scale (standard deviation). Normal distribution with mean = 0, standard deviation 1, has fraction <= -2, p = 0.0227501 Normal distribution with mean = 0, standard deviation 1, has fraction > -2, p = 0.97725 scale (standard deviation) = 0.647201 Normal distribution with mean = 0 has fraction <= -2, p = 0.001 Normal distribution with mean = 0 has fraction > -2, p = 0.999 Normal distribution with mean = 0.946339 has fraction <= -2 = 0.001 Normal distribution with mean = 0.946339 has fraction > -2 = 0.999