libs/math/doc/html/math_toolkit/stat_tut/weg/find_eg/find_scale_eg.html

   1 <html>
   2 <head>
   3 <meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
   4 <title>Find Scale (Standard Deviation) Example</title>
   5 <link rel="stylesheet" href="../../../../math.css" type="text/css">
   6 <meta name="generator" content="DocBook XSL Stylesheets V1.79.1">
   7 <link rel="home" href="../../../../index.html" title="Math Toolkit 2.11.0">
   8 <link rel="up" href="../find_eg.html" title="Find Location and Scale Examples">
   9 <link rel="prev" href="find_location_eg.html" title="Find Location (Mean) Example">
  10 <link rel="next" href="find_mean_and_sd_eg.html" title="Find mean and standard deviation example">
  11 </head>
  12 <body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
  13 <table cellpadding="2" width="100%"><tr>
  14 <td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../../../boost.png"></td>
  15 <td align="center"><a href="../../../../../../../../index.html">Home</a></td>
  16 <td align="center"><a href="../../../../../../../../libs/libraries.htm">Libraries</a></td>
  17 <td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
  18 <td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
  19 <td align="center"><a href="../../../../../../../../more/index.htm">More</a></td>
  20 </tr></table>
  21 <hr>
  22 <div class="spirit-nav">
  23 <a accesskey="p" href="find_location_eg.html"><img src="../../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../find_eg.html"><img src="../../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../../../index.html"><img src="../../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="find_mean_and_sd_eg.html"><img src="../../../../../../../../doc/src/images/next.png" alt="Next"></a>
  24 </div>
  25 <div class="section">
  26 <div class="titlepage"><div><div><h5 class="title">
  27 <a name="math_toolkit.stat_tut.weg.find_eg.find_scale_eg"></a><a class="link" href="find_scale_eg.html" title="Find Scale (Standard Deviation) Example">Find
  28           Scale (Standard Deviation) Example</a>
  29 </h5></div></div></div>
  30 <p>
  31             First we need some includes to access the <a class="link" href="../../../dist_ref/dists/normal_dist.html" title="Normal (Gaussian) Distribution">Normal
  32             Distribution</a>, the algorithms to find scale (and some std output
  33             of course).
  34           </p>
  35 <pre class="programlisting"><span class="preprocessor">#include</span> <span class="special">&lt;</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">distributions</span><span class="special">/</span><span class="identifier">normal</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">&gt;</span> <span class="comment">// for normal_distribution</span>
  36   <span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">normal</span><span class="special">;</span> <span class="comment">// typedef provides default type is double.</span>
  37 <span class="preprocessor">#include</span> <span class="special">&lt;</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">distributions</span><span class="special">/</span><span class="identifier">find_scale</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">&gt;</span>
  38   <span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">find_scale</span><span class="special">;</span>
  39   <span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">complement</span><span class="special">;</span> <span class="comment">// Needed if you want to use the complement version.</span>
  40   <span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">policies</span><span class="special">::</span><span class="identifier">policy</span><span class="special">;</span> <span class="comment">// Needed to specify the error handling policy.</span>
  41
  42 <span class="preprocessor">#include</span> <span class="special">&lt;</span><span class="identifier">iostream</span><span class="special">&gt;</span>
  43   <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">;</span> <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
  44 <span class="preprocessor">#include</span> <span class="special">&lt;</span><span class="identifier">iomanip</span><span class="special">&gt;</span>
  45   <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">setw</span><span class="special">;</span> <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">setprecision</span><span class="special">;</span>
  46 <span class="preprocessor">#include</span> <span class="special">&lt;</span><span class="identifier">limits</span><span class="special">&gt;</span>
  47   <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special">;</span>
  48 </pre>
  49 <p>
  50             For this example, we will use the standard <a class="link" href="../../../dist_ref/dists/normal_dist.html" title="Normal (Gaussian) Distribution">Normal
  51             Distribution</a>, with location (mean) zero and standard deviation
  52             (scale) unity. Conveniently, this is also the default for this implementation's
  53             constructor.
  54           </p>
  55 <pre class="programlisting"><span class="identifier">normal</span> <span class="identifier">N01</span><span class="special">;</span>  <span class="comment">// Default 'standard' normal distribution with zero mean</span>
  56 <span class="keyword">double</span> <span class="identifier">sd</span> <span class="special">=</span> <span class="number">1.</span><span class="special">;</span> <span class="comment">// and standard deviation is 1.</span>
  57 </pre>
  58 <p>
  59             Suppose we want to find a different normal distribution with standard
  60             deviation so that only fraction p (here 0.001 or 0.1%) are below a certain
  61             chosen limit (here -2. standard deviations).
  62           </p>
  63 <pre class="programlisting"><span class="keyword">double</span> <span class="identifier">z</span> <span class="special">=</span> <span class="special">-</span><span class="number">2.</span><span class="special">;</span> <span class="comment">// z to give prob p</span>
  64 <span class="keyword">double</span> <span class="identifier">p</span> <span class="special">=</span> <span class="number">0.001</span><span class="special">;</span> <span class="comment">// only 0.1% below z = -2</span>
  65
  66 <span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Normal distribution with mean = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">N01</span><span class="special">.</span><span class="identifier">location</span><span class="special">()</span>  <span class="comment">// aka N01.mean()</span>
  67   <span class="special">&lt;&lt;</span> <span class="string">", standard deviation "</span> <span class="special">&lt;&lt;</span> <span class="identifier">N01</span><span class="special">.</span><span class="identifier">scale</span><span class="special">()</span> <span class="comment">// aka N01.standard_deviation()</span>
  68   <span class="special">&lt;&lt;</span> <span class="string">", has "</span> <span class="special">&lt;&lt;</span> <span class="string">"fraction &lt;= "</span> <span class="special">&lt;&lt;</span> <span class="identifier">z</span>
  69   <span class="special">&lt;&lt;</span> <span class="string">", p = "</span>  <span class="special">&lt;&lt;</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">N01</span><span class="special">,</span> <span class="identifier">z</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
  70 <span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Normal distribution with mean = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">N01</span><span class="special">.</span><span class="identifier">location</span><span class="special">()</span>
  71   <span class="special">&lt;&lt;</span> <span class="string">", standard deviation "</span> <span class="special">&lt;&lt;</span> <span class="identifier">N01</span><span class="special">.</span><span class="identifier">scale</span><span class="special">()</span>
  72   <span class="special">&lt;&lt;</span> <span class="string">", has "</span> <span class="special">&lt;&lt;</span> <span class="string">"fraction &gt; "</span> <span class="special">&lt;&lt;</span> <span class="identifier">z</span>
  73   <span class="special">&lt;&lt;</span> <span class="string">", p = "</span>  <span class="special">&lt;&lt;</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">complement</span><span class="special">(</span><span class="identifier">N01</span><span class="special">,</span> <span class="identifier">z</span><span class="special">))</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// Note: uses complement.</span>
  74 </pre>
  75 <pre class="programlisting">Normal distribution with mean = 0 has fraction &lt;= -2, p = 0.0227501
  76 Normal distribution with mean = 0 has fraction &gt; -2, p = 0.97725
  77 </pre>
  78 <p>
  79             Noting that p = 0.02 instead of our target of 0.001, we can now use
  80             <code class="computeroutput"><span class="identifier">find_scale</span></code> to give a
  81             new standard deviation.
  82           </p>
  83 <pre class="programlisting"><span class="keyword">double</span> <span class="identifier">l</span> <span class="special">=</span> <span class="identifier">N01</span><span class="special">.</span><span class="identifier">location</span><span class="special">();</span>
  84 <span class="keyword">double</span> <span class="identifier">s</span> <span class="special">=</span> <span class="identifier">find_scale</span><span class="special">&lt;</span><span class="identifier">normal</span><span class="special">&gt;(</span><span class="identifier">z</span><span class="special">,</span> <span class="identifier">p</span><span class="special">,</span> <span class="identifier">l</span><span class="special">);</span>
  85 <span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"scale (standard deviation) = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">s</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
  86 </pre>
  87 <p>
  88             that outputs:
  89           </p>
  90 <pre class="programlisting">scale (standard deviation) = 0.647201
  91 </pre>
  92 <p>
  93             showing that we need to reduce the standard deviation from 1. to 0.65.
  94           </p>
  95 <p>
  96             Then we can check that we have achieved our objective by constructing
  97             a new distribution with the new standard deviation (but same zero mean):
  98           </p>
  99 <pre class="programlisting"><span class="identifier">normal</span> <span class="identifier">np001pc</span><span class="special">(</span><span class="identifier">N01</span><span class="special">.</span><span class="identifier">location</span><span class="special">(),</span> <span class="identifier">s</span><span class="special">);</span>
 100 </pre>
 101 <p>
 102             And re-calculating the fraction below (and above) our chosen limit.
 103           </p>
 104 <pre class="programlisting"><span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Normal distribution with mean = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">l</span>
 105   <span class="special">&lt;&lt;</span> <span class="string">" has "</span> <span class="special">&lt;&lt;</span> <span class="string">"fraction &lt;= "</span> <span class="special">&lt;&lt;</span> <span class="identifier">z</span>
 106   <span class="special">&lt;&lt;</span> <span class="string">", p = "</span>  <span class="special">&lt;&lt;</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">np001pc</span><span class="special">,</span> <span class="identifier">z</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
 107 <span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Normal distribution with mean = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">l</span>
 108   <span class="special">&lt;&lt;</span> <span class="string">" has "</span> <span class="special">&lt;&lt;</span> <span class="string">"fraction &gt; "</span> <span class="special">&lt;&lt;</span> <span class="identifier">z</span>
 109   <span class="special">&lt;&lt;</span> <span class="string">", p = "</span>  <span class="special">&lt;&lt;</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">complement</span><span class="special">(</span><span class="identifier">np001pc</span><span class="special">,</span> <span class="identifier">z</span><span class="special">))</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
 110 </pre>
 111 <pre class="programlisting">Normal distribution with mean = 0 has fraction &lt;= -2, p = 0.001
 112 Normal distribution with mean = 0 has fraction &gt; -2, p = 0.999
 113 </pre>
 114 <h5>
 115 <a name="math_toolkit.stat_tut.weg.find_eg.find_scale_eg.h0"></a>
 116             <span class="phrase"><a name="math_toolkit.stat_tut.weg.find_eg.find_scale_eg.controlling_how_errors_from_find"></a></span><a class="link" href="find_scale_eg.html#math_toolkit.stat_tut.weg.find_eg.find_scale_eg.controlling_how_errors_from_find">Controlling
 117             how Errors from find_scale are handled</a>
 118           </h5>
 119 <p>
 120             We can also control the policy for handling various errors. For example,
 121             we can define a new (possibly unwise) policy to ignore domain errors
 122             ('bad' arguments).
 123           </p>
 124 <p>
 125             Unless we are using the boost::math namespace, we will need:
 126           </p>
 127 <pre class="programlisting"><span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">policies</span><span class="special">::</span><span class="identifier">policy</span><span class="special">;</span>
 128 <span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">policies</span><span class="special">::</span><span class="identifier">domain_error</span><span class="special">;</span>
 129 <span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">policies</span><span class="special">::</span><span class="identifier">ignore_error</span><span class="special">;</span>
 130 </pre>
 131 <p>
 132             Using a typedef is convenient, especially if it is re-used, although
 133             it is not required, as the various examples below show.
 134           </p>
 135 <pre class="programlisting"><span class="keyword">typedef</span> <span class="identifier">policy</span><span class="special">&lt;</span><span class="identifier">domain_error</span><span class="special">&lt;</span><span class="identifier">ignore_error</span><span class="special">&gt;</span> <span class="special">&gt;</span> <span class="identifier">ignore_domain_policy</span><span class="special">;</span>
 136 <span class="comment">// find_scale with new policy, using typedef.</span>
 137 <span class="identifier">l</span> <span class="special">=</span> <span class="identifier">find_scale</span><span class="special">&lt;</span><span class="identifier">normal</span><span class="special">&gt;(</span><span class="identifier">z</span><span class="special">,</span> <span class="identifier">p</span><span class="special">,</span> <span class="identifier">l</span><span class="special">,</span> <span class="identifier">ignore_domain_policy</span><span class="special">());</span>
 138 <span class="comment">// Default policy policy&lt;&gt;, needs using boost::math::policies::policy;</span>
 139
 140 <span class="identifier">l</span> <span class="special">=</span> <span class="identifier">find_scale</span><span class="special">&lt;</span><span class="identifier">normal</span><span class="special">&gt;(</span><span class="identifier">z</span><span class="special">,</span> <span class="identifier">p</span><span class="special">,</span> <span class="identifier">l</span><span class="special">,</span> <span class="identifier">policy</span><span class="special">&lt;&gt;());</span>
 141 <span class="comment">// Default policy, fully specified.</span>
 142 <span class="identifier">l</span> <span class="special">=</span> <span class="identifier">find_scale</span><span class="special">&lt;</span><span class="identifier">normal</span><span class="special">&gt;(</span><span class="identifier">z</span><span class="special">,</span> <span class="identifier">p</span><span class="special">,</span> <span class="identifier">l</span><span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">policies</span><span class="special">::</span><span class="identifier">policy</span><span class="special">&lt;&gt;());</span>
 143 <span class="comment">// New policy, without typedef.</span>
 144 <span class="identifier">l</span> <span class="special">=</span> <span class="identifier">find_scale</span><span class="special">&lt;</span><span class="identifier">normal</span><span class="special">&gt;(</span><span class="identifier">z</span><span class="special">,</span> <span class="identifier">p</span><span class="special">,</span> <span class="identifier">l</span><span class="special">,</span> <span class="identifier">policy</span><span class="special">&lt;</span><span class="identifier">domain_error</span><span class="special">&lt;</span><span class="identifier">ignore_error</span><span class="special">&gt;</span> <span class="special">&gt;());</span>
 145 </pre>
 146 <p>
 147             If we want to express a probability, say 0.999, that is a complement,
 148             <code class="computeroutput"><span class="number">1</span> <span class="special">-</span>
 149             <span class="identifier">p</span></code> we should not even think
 150             of writing <code class="computeroutput"><span class="identifier">find_scale</span><span class="special">&lt;</span><span class="identifier">normal</span><span class="special">&gt;(</span><span class="identifier">z</span><span class="special">,</span> <span class="number">1</span> <span class="special">-</span>
 151             <span class="identifier">p</span><span class="special">,</span>
 152             <span class="identifier">l</span><span class="special">)</span></code>,
 153             but use the <a class="link" href="../../overview/complements.html" title="Complements are supported too - and when to use them">complements</a>
 154             version (see <a class="link" href="../../overview/complements.html#why_complements">why complements?</a>).
 155           </p>
 156 <pre class="programlisting"><span class="identifier">z</span> <span class="special">=</span> <span class="special">-</span><span class="number">2.</span><span class="special">;</span>
 157 <span class="keyword">double</span> <span class="identifier">q</span> <span class="special">=</span> <span class="number">0.999</span><span class="special">;</span> <span class="comment">// = 1 - p; // complement of 0.001.</span>
 158 <span class="identifier">sd</span> <span class="special">=</span> <span class="identifier">find_scale</span><span class="special">&lt;</span><span class="identifier">normal</span><span class="special">&gt;(</span><span class="identifier">complement</span><span class="special">(</span><span class="identifier">z</span><span class="special">,</span> <span class="identifier">q</span><span class="special">,</span> <span class="identifier">l</span><span class="special">));</span>
 159
 160 <span class="identifier">normal</span> <span class="identifier">np95pc</span><span class="special">(</span><span class="identifier">l</span><span class="special">,</span> <span class="identifier">sd</span><span class="special">);</span> <span class="comment">// Same standard_deviation (scale) but with mean(scale) shifted</span>
 161 <span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Normal distribution with mean = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">l</span> <span class="special">&lt;&lt;</span> <span class="string">" has "</span>
 162   <span class="special">&lt;&lt;</span> <span class="string">"fraction &lt;= "</span> <span class="special">&lt;&lt;</span> <span class="identifier">z</span> <span class="special">&lt;&lt;</span> <span class="string">" = "</span>  <span class="special">&lt;&lt;</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">np95pc</span><span class="special">,</span> <span class="identifier">z</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
 163 <span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Normal distribution with mean = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">l</span> <span class="special">&lt;&lt;</span> <span class="string">" has "</span>
 164   <span class="special">&lt;&lt;</span> <span class="string">"fraction &gt; "</span> <span class="special">&lt;&lt;</span> <span class="identifier">z</span> <span class="special">&lt;&lt;</span> <span class="string">" = "</span>  <span class="special">&lt;&lt;</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">complement</span><span class="special">(</span><span class="identifier">np95pc</span><span class="special">,</span> <span class="identifier">z</span><span class="special">))</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
 165 </pre>
 166 <p>
 167             Sadly, it is all too easy to get probabilities the wrong way round, when
 168             you may get a warning like this:
 169           </p>
 170 <pre class="programlisting">Message from thrown exception was:
 171    Error in function boost::math::find_scale&lt;Dist, Policy&gt;(complement(double, double, double, Policy)):
 172    Computed scale (-0.48043523852179076) is &lt;= 0! Was the complement intended?
 173 </pre>
 174 <p>
 175             The default error handling policy is to throw an exception with this
 176             message, but if you chose a policy to ignore the error, the (impossible)
 177             negative scale is quietly returned.
 178           </p>
 179 <p>
 180             See <a href="../../../../../../example/find_scale_example.cpp" target="_top">find_scale_example.cpp</a>
 181             for full source code: the program output looks like this:
 182           </p>
 183 <pre class="programlisting">Example: Find scale (standard deviation).
 184 Normal distribution with mean = 0, standard deviation 1, has fraction &lt;= -2, p = 0.0227501
 185 Normal distribution with mean = 0, standard deviation 1, has fraction &gt; -2, p = 0.97725
 186 scale (standard deviation) = 0.647201
 187 Normal distribution with mean = 0 has fraction &lt;= -2, p = 0.001
 188 Normal distribution with mean = 0 has fraction &gt; -2, p = 0.999
 189 Normal distribution with mean = 0.946339 has fraction &lt;= -2 = 0.001
 190 Normal distribution with mean = 0.946339 has fraction &gt; -2 = 0.999
 191 </pre>
 192 </div>
 193 <table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
 194 <td align="left"></td>
 195 <td align="right"><div class="copyright-footer">Copyright &#169; 2006-2019 Nikhar
 196       Agrawal, Anton Bikineev, Paul A. Bristow, Marco Guazzone, Christopher Kormanyos,
 197       Hubert Holin, Bruno Lalande, John Maddock, Jeremy Murphy, Matthew Pulver, Johan
 198       R&#229;de, Gautam Sewani, Benjamin Sobotta, Nicholas Thompson, Thijs van den Berg,
 199       Daryle Walker and Xiaogang Zhang<p>
 200         Distributed under the Boost Software License, Version 1.0. (See accompanying
 201         file LICENSE_1_0.txt or copy at <a href="http://www.boost.org/LICENSE_1_0.txt" target="_top">http://www.boost.org/LICENSE_1_0.txt</a>)
 202       </p>
 203 </div></td>
 204 </tr></table>
 205 <hr>
 206 <div class="spirit-nav">
 207 <a accesskey="p" href="find_location_eg.html"><img src="../../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../find_eg.html"><img src="../../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../../../index.html"><img src="../../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="find_mean_and_sd_eg.html"><img src="../../../../../../../../doc/src/images/next.png" alt="Next"></a>
 208 </div>
 209 </body>
 210 </html>