Find Scale (Standard Deviation) Example

| Home | Libraries | People | FAQ | More

###### Find Scale (Standard Deviation) Example

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

##### Controlling how Errors from find_scale are handled

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