blob: 16797c15b134db23b6b23385003a4906faf03971 [file] [log] [blame]
<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
<title>Some Miscellaneous Examples of the Normal (Gaussian) Distribution</title>
<link rel="stylesheet" href="../../../../../../../../../../doc/src/boostbook.css" type="text/css">
<meta name="generator" content="DocBook XSL Stylesheets V1.74.0">
<link rel="home" href="../../../../../index.html" title="Math Toolkit">
<link rel="up" href="../normal_example.html" title="Normal Distribution Examples">
<link rel="prev" href="../normal_example.html" title="Normal Distribution Examples">
<link rel="next" href="../nccs_eg.html" title="Non Central Chi Squared Example">
</head>
<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
<table cellpadding="2" width="100%"><tr>
<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../../../../../boost.png"></td>
<td align="center"><a href="../../../../../../../../../../index.html">Home</a></td>
<td align="center"><a href="../../../../../../../../../../libs/libraries.htm">Libraries</a></td>
<td align="center"><a href="http://www.boost.org/users/people.html">People</a></td>
<td align="center"><a href="http://www.boost.org/users/faq.html">FAQ</a></td>
<td align="center"><a href="../../../../../../../../../../more/index.htm">More</a></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="../normal_example.html"><img src="../../../../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../normal_example.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="../nccs_eg.html"><img src="../../../../../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
<div class="section" lang="en">
<div class="titlepage"><div><div><h6 class="title">
<a name="math_toolkit.dist.stat_tut.weg.normal_example.normal_misc"></a><a class="link" href="normal_misc.html" title="Some Miscellaneous Examples of the Normal (Gaussian) Distribution">
Some Miscellaneous Examples of the Normal (Gaussian) Distribution</a>
</h6></div></div></div>
<p>
The sample program <a href="../../../../../../../../example/normal_misc_examples.cpp" target="_top">normal_misc_examples.cpp</a>
illustrates their use.
</p>
<a name="math_toolkit.dist.stat_tut.weg.normal_example.normal_misc.traditional_tables"></a><h5>
<a name="id972622"></a>
<a class="link" href="normal_misc.html#math_toolkit.dist.stat_tut.weg.normal_example.normal_misc.traditional_tables">Traditional
Tables</a>
</h5>
<p>
</p>
<p>
First we need some includes to access the normal distribution (and
some std output of course).
</p>
<p>
</p>
<p>
</p>
<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> <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>
<span class="preprocessor">#include</span> <span class="special">&lt;</span><span class="identifier">iostream</span><span class="special">&gt;</span>
<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> <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">left</span><span class="special">;</span> <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">showpoint</span><span class="special">;</span> <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">noshowpoint</span><span class="special">;</span>
<span class="preprocessor">#include</span> <span class="special">&lt;</span><span class="identifier">iomanip</span><span class="special">&gt;</span>
<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>
<span class="preprocessor">#include</span> <span class="special">&lt;</span><span class="identifier">limits</span><span class="special">&gt;</span>
<span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special">;</span>
<span class="keyword">int</span> <span class="identifier">main</span><span class="special">()</span>
<span class="special">{</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Example: Normal distribution, Miscellaneous Applications."</span><span class="special">;</span>
<span class="keyword">try</span>
<span class="special">{</span>
<span class="special">{</span> <span class="comment">// Traditional tables and values.</span></pre>
<p>
</p>
<p>
</p>
<p>
Let's start by printing some traditional tables.
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">step</span> <span class="special">=</span> <span class="number">1.</span><span class="special">;</span> <span class="comment">// in z
</span><span class="keyword">double</span> <span class="identifier">range</span> <span class="special">=</span> <span class="number">4</span><span class="special">;</span> <span class="comment">// min and max z = -range to +range.
</span><span class="keyword">int</span> <span class="identifier">precision</span> <span class="special">=</span> <span class="number">17</span><span class="special">;</span> <span class="comment">// traditional tables are only computed to much lower precision.
</span>
<span class="comment">// Construct a standard normal distribution s
</span> <span class="identifier">normal</span> <span class="identifier">s</span><span class="special">;</span> <span class="comment">// (default mean = zero, and standard deviation = unity)
</span> <span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Standard normal distribution, mean = "</span><span class="special">&lt;&lt;</span> <span class="identifier">s</span><span class="special">.</span><span class="identifier">mean</span><span class="special">()</span>
<span class="special">&lt;&lt;</span> <span class="string">", standard deviation = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">s</span><span class="special">.</span><span class="identifier">standard_deviation</span><span class="special">()</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span></pre>
<p>
</p>
<p>
</p>
<p>
First the probability distribution function (pdf).
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"><span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Probability distribution function values"</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">" z "</span> <span class="string">" pdf "</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="number">5</span><span class="special">);</span>
<span class="keyword">for</span> <span class="special">(</span><span class="keyword">double</span> <span class="identifier">z</span> <span class="special">=</span> <span class="special">-</span><span class="identifier">range</span><span class="special">;</span> <span class="identifier">z</span> <span class="special">&lt;</span> <span class="identifier">range</span> <span class="special">+</span> <span class="identifier">step</span><span class="special">;</span> <span class="identifier">z</span> <span class="special">+=</span> <span class="identifier">step</span><span class="special">)</span>
<span class="special">{</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="identifier">left</span> <span class="special">&lt;&lt;</span> <span class="identifier">setprecision</span><span class="special">(</span><span class="number">3</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="identifier">setw</span><span class="special">(</span><span class="number">6</span><span class="special">)</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">setprecision</span><span class="special">(</span><span class="identifier">precision</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="identifier">setw</span><span class="special">(</span><span class="number">12</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="identifier">pdf</span><span class="special">(</span><span class="identifier">s</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="special">}</span>
<span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="number">6</span><span class="special">);</span> <span class="comment">// default</span></pre>
<p>
</p>
<p>
</p>
<p>
And the area under the normal curve from -&#8734; up to z, the cumulative
distribution function (cdf).
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"><span class="comment">// For a standard normal distribution
</span><span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Standard normal mean = "</span><span class="special">&lt;&lt;</span> <span class="identifier">s</span><span class="special">.</span><span class="identifier">mean</span><span class="special">()</span>
<span class="special">&lt;&lt;</span> <span class="string">", standard deviation = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">s</span><span class="special">.</span><span class="identifier">standard_deviation</span><span class="special">()</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Integral (area under the curve) from - infinity up to z "</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">" z "</span> <span class="string">" cdf "</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="keyword">for</span> <span class="special">(</span><span class="keyword">double</span> <span class="identifier">z</span> <span class="special">=</span> <span class="special">-</span><span class="identifier">range</span><span class="special">;</span> <span class="identifier">z</span> <span class="special">&lt;</span> <span class="identifier">range</span> <span class="special">+</span> <span class="identifier">step</span><span class="special">;</span> <span class="identifier">z</span> <span class="special">+=</span> <span class="identifier">step</span><span class="special">)</span>
<span class="special">{</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="identifier">left</span> <span class="special">&lt;&lt;</span> <span class="identifier">setprecision</span><span class="special">(</span><span class="number">3</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="identifier">setw</span><span class="special">(</span><span class="number">6</span><span class="special">)</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">setprecision</span><span class="special">(</span><span class="identifier">precision</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="identifier">setw</span><span class="special">(</span><span class="number">12</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">s</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="special">}</span>
<span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="number">6</span><span class="special">);</span> <span class="comment">// default</span></pre>
<p>
</p>
<p>
</p>
<p>
And all this you can do with a nanoscopic amount of work compared
to the team of <span class="bold"><strong>human computers</strong></span> toiling
with Milton Abramovitz and Irene Stegen at the US National Bureau
of Standards (now <a href="http://www.nist.gov" target="_top">NIST</a>).
Starting in 1938, their "Handbook of Mathematical Functions
with Formulas, Graphs and Mathematical Tables", was eventually
published in 1964, and has been reprinted numerous times since. (A
major replacement is planned at <a href="http://dlmf.nist.gov" target="_top">Digital
Library of Mathematical Functions</a>).
</p>
<p>
</p>
<p>
Pretty-printing a traditional 2-dimensional table is left as an exercise
for the student, but why bother now that the Math Toolkit lets you
write
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">z</span> <span class="special">=</span> <span class="number">2.</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Area for z = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">z</span> <span class="special">&lt;&lt;</span> <span class="string">" is "</span> <span class="special">&lt;&lt;</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">s</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">// to get the area for z.</span></pre>
<p>
</p>
<p>
</p>
<p>
Correspondingly, we can obtain the traditional 'critical' values
for significance levels. For the 95% confidence level, the significance
level usually called alpha, is 0.05 = 1 - 0.95 (for a one-sided test),
so we can write
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"> <span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"95% of area has a z below "</span> <span class="special">&lt;&lt;</span> <span class="identifier">quantile</span><span class="special">(</span><span class="identifier">s</span><span class="special">,</span> <span class="number">0.95</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="comment">// 95% of area has a z below 1.64485</span></pre>
<p>
</p>
<p>
</p>
<p>
and a two-sided test (a comparison between two levels, rather than
a one-sided test)
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"> <span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"95% of area has a z between "</span> <span class="special">&lt;&lt;</span> <span class="identifier">quantile</span><span class="special">(</span><span class="identifier">s</span><span class="special">,</span> <span class="number">0.975</span><span class="special">)</span>
<span class="special">&lt;&lt;</span> <span class="string">" and "</span> <span class="special">&lt;&lt;</span> <span class="special">-</span><span class="identifier">quantile</span><span class="special">(</span><span class="identifier">s</span><span class="special">,</span> <span class="number">0.975</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="comment">// 95% of area has a z between 1.95996 and -1.95996</span></pre>
<p>
</p>
<p>
</p>
<p>
First, define a table of significance levels: these are the probabilities
that the true occurrence frequency lies outside the calculated interval.
</p>
<p>
</p>
<p>
It is convenient to have an alpha level for the probability that
z lies outside just one standard deviation. This will not be some
nice neat number like 0.05, but we can easily calculate it,
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">alpha1</span> <span class="special">=</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">s</span><span class="special">,</span> <span class="special">-</span><span class="number">1</span><span class="special">)</span> <span class="special">*</span> <span class="number">2</span><span class="special">;</span> <span class="comment">// 0.3173105078629142
</span><span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="identifier">setprecision</span><span class="special">(</span><span class="number">17</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="string">"Significance level for z == 1 is "</span> <span class="special">&lt;&lt;</span> <span class="identifier">alpha1</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span></pre>
<p>
</p>
<p>
</p>
<p>
and place in our array of favorite alpha values.
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">alpha</span><span class="special">[]</span> <span class="special">=</span> <span class="special">{</span><span class="number">0.3173105078629142</span><span class="special">,</span> <span class="comment">// z for 1 standard deviation.
</span> <span class="number">0.20</span><span class="special">,</span> <span class="number">0.1</span><span class="special">,</span> <span class="number">0.05</span><span class="special">,</span> <span class="number">0.01</span><span class="special">,</span> <span class="number">0.001</span><span class="special">,</span> <span class="number">0.0001</span><span class="special">,</span> <span class="number">0.00001</span> <span class="special">};</span></pre>
<p>
</p>
<p>
</p>
<p>
Confidence value as % is (1 - alpha) * 100 (so alpha 0.05 == 95%
confidence) that the true occurrence frequency lies <span class="bold"><strong>inside</strong></span>
the calculated interval.
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"><span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"level of significance (alpha)"</span> <span class="special">&lt;&lt;</span> <span class="identifier">setprecision</span><span class="special">(</span><span class="number">4</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"2-sided 1 -sided z(alpha) "</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="keyword">for</span> <span class="special">(</span><span class="keyword">int</span> <span class="identifier">i</span> <span class="special">=</span> <span class="number">0</span><span class="special">;</span> <span class="identifier">i</span> <span class="special">&lt;</span> <span class="keyword">sizeof</span><span class="special">(</span><span class="identifier">alpha</span><span class="special">)/</span><span class="keyword">sizeof</span><span class="special">(</span><span class="identifier">alpha</span><span class="special">[</span><span class="number">0</span><span class="special">]);</span> <span class="special">++</span><span class="identifier">i</span><span class="special">)</span>
<span class="special">{</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="identifier">setw</span><span class="special">(</span><span class="number">15</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="identifier">alpha</span><span class="special">[</span><span class="identifier">i</span><span class="special">]</span> <span class="special">&lt;&lt;</span> <span class="identifier">setw</span><span class="special">(</span><span class="number">15</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="identifier">alpha</span><span class="special">[</span><span class="identifier">i</span><span class="special">]</span> <span class="special">/</span><span class="number">2</span> <span class="special">&lt;&lt;</span> <span class="identifier">setw</span><span class="special">(</span><span class="number">10</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="identifier">quantile</span><span class="special">(</span><span class="identifier">complement</span><span class="special">(</span><span class="identifier">s</span><span class="special">,</span> <span class="identifier">alpha</span><span class="special">[</span><span class="identifier">i</span><span class="special">]/</span><span class="number">2</span><span class="special">))</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="comment">// Use quantile(complement(s, alpha[i]/2)) to avoid potential loss of accuracy from quantile(s, 1 - alpha[i]/2)
</span><span class="special">}</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span></pre>
<p>
</p>
<p>
</p>
<p>
Notice the distinction between one-sided (also called one-tailed)
where we are using a &gt; <span class="bold"><strong>or</strong></span> &lt;
test (and not both) and considering the area of the tail (integral)
from z up to +&#8734;, and a two-sided test where we are using two &gt;
<span class="bold"><strong>and</strong></span> &lt; tests, and thus considering
two tails, from -&#8734; up to z low and z high up to +&#8734;.
</p>
<p>
</p>
<p>
So the 2-sided values alpha[i] are calculated using alpha[i]/2.
</p>
<p>
</p>
<p>
If we consider a simple example of alpha = 0.05, then for a two-sided
test, the lower tail area from -&#8734; up to -1.96 is 0.025 (alpha/2) and
the upper tail area from +z up to +1.96 is also 0.025 (alpha/2),
and the area between -1.96 up to 12.96 is alpha = 0.95. and the sum
of the two tails is 0.025 + 0.025 = 0.05,
</p>
<p>
</p>
<a name="math_toolkit.dist.stat_tut.weg.normal_example.normal_misc.standard_deviations_either_side_of_the_mean"></a><h5>
<a name="id976060"></a>
<a class="link" href="normal_misc.html#math_toolkit.dist.stat_tut.weg.normal_example.normal_misc.standard_deviations_either_side_of_the_mean">Standard
deviations either side of the Mean</a>
</h5>
<p>
</p>
<p>
Armed with the cumulative distribution function, we can easily calculate
the easy to remember proportion of values that lie within 1, 2 and
3 standard deviations from the mean.
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="number">3</span><span class="special">);</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="identifier">showpoint</span> <span class="special">&lt;&lt;</span> <span class="string">"cdf(s, s.standard_deviation()) = "</span>
<span class="special">&lt;&lt;</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">s</span><span class="special">,</span> <span class="identifier">s</span><span class="special">.</span><span class="identifier">standard_deviation</span><span class="special">())</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// from -infinity to 1 sd
</span><span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"cdf(complement(s, s.standard_deviation())) = "</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">s</span><span class="special">,</span> <span class="identifier">s</span><span class="special">.</span><span class="identifier">standard_deviation</span><span class="special">()))</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Fraction 1 standard deviation within either side of mean is "</span>
<span class="special">&lt;&lt;</span> <span class="number">1</span> <span class="special">-</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">complement</span><span class="special">(</span><span class="identifier">s</span><span class="special">,</span> <span class="identifier">s</span><span class="special">.</span><span class="identifier">standard_deviation</span><span class="special">()))</span> <span class="special">*</span> <span class="number">2</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Fraction 2 standard deviations within either side of mean is "</span>
<span class="special">&lt;&lt;</span> <span class="number">1</span> <span class="special">-</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">complement</span><span class="special">(</span><span class="identifier">s</span><span class="special">,</span> <span class="number">2</span> <span class="special">*</span> <span class="identifier">s</span><span class="special">.</span><span class="identifier">standard_deviation</span><span class="special">()))</span> <span class="special">*</span> <span class="number">2</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Fraction 3 standard deviations within either side of mean is "</span>
<span class="special">&lt;&lt;</span> <span class="number">1</span> <span class="special">-</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">complement</span><span class="special">(</span><span class="identifier">s</span><span class="special">,</span> <span class="number">3</span> <span class="special">*</span> <span class="identifier">s</span><span class="special">.</span><span class="identifier">standard_deviation</span><span class="special">()))</span> <span class="special">*</span> <span class="number">2</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span></pre>
<p>
</p>
<p>
</p>
<p>
To a useful precision, the 1, 2 &amp; 3 percentages are 68, 95 and
99.7, and these are worth memorising as useful 'rules of thumb',
as, for example, in <a href="http://en.wikipedia.org/wiki/Standard_deviation" target="_top">standard
deviation</a>:
</p>
<p>
</p>
<pre class="programlisting">Fraction 1 standard deviation within either side of mean is 0.683
Fraction 2 standard deviations within either side of mean is 0.954
Fraction 3 standard deviations within either side of mean is 0.997
</pre>
<p>
</p>
<p>
We could of course get some really accurate values for these <a href="http://en.wikipedia.org/wiki/Confidence_interval" target="_top">confidence
intervals</a> by using cout.precision(15);
</p>
<p>
</p>
<pre class="programlisting">Fraction 1 standard deviation within either side of mean is 0.682689492137086
Fraction 2 standard deviations within either side of mean is 0.954499736103642
Fraction 3 standard deviations within either side of mean is 0.997300203936740
</pre>
<p>
</p>
<p>
But before you get too excited about this impressive precision, don't
forget that the <span class="bold"><strong>confidence intervals of the
standard deviation</strong></span> are surprisingly wide, especially if
you have estimated the standard deviation from only a few measurements.
</p>
<p>
</p>
<a name="math_toolkit.dist.stat_tut.weg.normal_example.normal_misc.some_simple_examples"></a><h5>
<a name="id976885"></a>
<a class="link" href="normal_misc.html#math_toolkit.dist.stat_tut.weg.normal_example.normal_misc.some_simple_examples">Some
simple examples</a>
</h5>
<a name="math_toolkit.dist.stat_tut.weg.normal_example.normal_misc.life_of_light_bulbs"></a><h5>
<a name="id976898"></a>
<a class="link" href="normal_misc.html#math_toolkit.dist.stat_tut.weg.normal_example.normal_misc.life_of_light_bulbs">Life
of light bulbs</a>
</h5>
<p>
</p>
<p>
Examples from K. Krishnamoorthy, Handbook of Statistical Distributions
with Applications, ISBN 1 58488 635 8, page 125... implemented using
the Math Toolkit library.
</p>
<p>
</p>
<p>
A few very simple examples are shown here:
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"><span class="comment">// K. Krishnamoorthy, Handbook of Statistical Distributions with Applications,
</span> <span class="comment">// ISBN 1 58488 635 8, page 125, example 10.3.5</span></pre>
<p>
</p>
<p>
</p>
<p>
Mean lifespan of 100 W bulbs is 1100 h with standard deviation of
100 h. Assuming, perhaps with little evidence and much faith, that
the distribution is normal, we construct a normal distribution called
<span class="emphasis"><em>bulbs</em></span> with these values:
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">mean_life</span> <span class="special">=</span> <span class="number">1100.</span><span class="special">;</span>
<span class="keyword">double</span> <span class="identifier">life_standard_deviation</span> <span class="special">=</span> <span class="number">100.</span><span class="special">;</span>
<span class="identifier">normal</span> <span class="identifier">bulbs</span><span class="special">(</span><span class="identifier">mean_life</span><span class="special">,</span> <span class="identifier">life_standard_deviation</span><span class="special">);</span>
<span class="keyword">double</span> <span class="identifier">expected_life</span> <span class="special">=</span> <span class="number">1000.</span><span class="special">;</span></pre>
<p>
</p>
<p>
</p>
<p>
The we can use the Cumulative distribution function to predict fractions
(or percentages, if * 100) that will last various lifetimes.
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"><span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Fraction of bulbs that will last at best (&lt;=) "</span> <span class="comment">// P(X &lt;= 1000)
</span> <span class="special">&lt;&lt;</span> <span class="identifier">expected_life</span> <span class="special">&lt;&lt;</span> <span class="string">" is "</span><span class="special">&lt;&lt;</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">bulbs</span><span class="special">,</span> <span class="identifier">expected_life</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Fraction of bulbs that will last at least (&gt;) "</span> <span class="comment">// P(X &gt; 1000)
</span> <span class="special">&lt;&lt;</span> <span class="identifier">expected_life</span> <span class="special">&lt;&lt;</span> <span class="string">" is "</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">bulbs</span><span class="special">,</span> <span class="identifier">expected_life</span><span class="special">))</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="keyword">double</span> <span class="identifier">min_life</span> <span class="special">=</span> <span class="number">900</span><span class="special">;</span>
<span class="keyword">double</span> <span class="identifier">max_life</span> <span class="special">=</span> <span class="number">1200</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Fraction of bulbs that will last between "</span>
<span class="special">&lt;&lt;</span> <span class="identifier">min_life</span> <span class="special">&lt;&lt;</span> <span class="string">" and "</span> <span class="special">&lt;&lt;</span> <span class="identifier">max_life</span> <span class="special">&lt;&lt;</span> <span class="string">" is "</span>
<span class="special">&lt;&lt;</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">bulbs</span><span class="special">,</span> <span class="identifier">max_life</span><span class="special">)</span> <span class="comment">// P(X &lt;= 1200)
</span> <span class="special">-</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">bulbs</span><span class="special">,</span> <span class="identifier">min_life</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// P(X &lt;= 900)</span></pre>
<p>
</p>
<p>
</p>
<div class="note"><table border="0" summary="Note">
<tr>
<td rowspan="2" align="center" valign="top" width="25"><img alt="[Note]" src="../../../../../../../../../../doc/src/images/note.png"></td>
<th align="left">Note</th>
</tr>
<tr><td align="left" valign="top"><p>
Real-life failures are often very ab-normal, with a significant
number that 'dead-on-arrival' or suffer failure very early in their
life: the lifetime of the survivors of 'early mortality' may be
well described by the normal distribution.
</p></td></tr>
</table></div>
<p>
</p>
<a name="math_toolkit.dist.stat_tut.weg.normal_example.normal_misc.how_many_onions_"></a><h5>
<a name="id977409"></a>
<a class="link" href="normal_misc.html#math_toolkit.dist.stat_tut.weg.normal_example.normal_misc.how_many_onions_">How
many onions?</a>
</h5>
<p>
</p>
<p>
Weekly demand for 5 lb sacks of onions at a store is normally distributed
with mean 140 sacks and standard deviation 10.
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">mean</span> <span class="special">=</span> <span class="number">140.</span><span class="special">;</span> <span class="comment">// sacks per week.
</span><span class="keyword">double</span> <span class="identifier">standard_deviation</span> <span class="special">=</span> <span class="number">10</span><span class="special">;</span>
<span class="identifier">normal</span> <span class="identifier">sacks</span><span class="special">(</span><span class="identifier">mean</span><span class="special">,</span> <span class="identifier">standard_deviation</span><span class="special">);</span>
<span class="keyword">double</span> <span class="identifier">stock</span> <span class="special">=</span> <span class="number">160.</span><span class="special">;</span> <span class="comment">// per week.
</span><span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Percentage of weeks overstocked "</span>
<span class="special">&lt;&lt;</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">sacks</span><span class="special">,</span> <span class="identifier">stock</span><span class="special">)</span> <span class="special">*</span> <span class="number">100.</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// P(X &lt;=160)
</span><span class="comment">// Percentage of weeks overstocked 97.7</span></pre>
<p>
</p>
<p>
</p>
<p>
So there will be lots of mouldy onions! So we should be able to say
what stock level will meet demand 95% of the weeks.
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">stock_95</span> <span class="special">=</span> <span class="identifier">quantile</span><span class="special">(</span><span class="identifier">sacks</span><span class="special">,</span> <span class="number">0.95</span><span class="special">);</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Store should stock "</span> <span class="special">&lt;&lt;</span> <span class="keyword">int</span><span class="special">(</span><span class="identifier">stock_95</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="string">" sacks to meet 95% of demands."</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span></pre>
<p>
</p>
<p>
</p>
<p>
And it is easy to estimate how to meet 80% of demand, and waste even
less.
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">stock_80</span> <span class="special">=</span> <span class="identifier">quantile</span><span class="special">(</span><span class="identifier">sacks</span><span class="special">,</span> <span class="number">0.80</span><span class="special">);</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Store should stock "</span> <span class="special">&lt;&lt;</span> <span class="keyword">int</span><span class="special">(</span><span class="identifier">stock_80</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="string">" sacks to meet 8 out of 10 demands."</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
</pre>
<p>
</p>
<p>
</p>
<a name="math_toolkit.dist.stat_tut.weg.normal_example.normal_misc.packing_beef"></a><h5>
<a name="id977827"></a>
<a class="link" href="normal_misc.html#math_toolkit.dist.stat_tut.weg.normal_example.normal_misc.packing_beef">Packing
beef</a>
</h5>
<p>
</p>
<p>
A machine is set to pack 3 kg of ground beef per pack. Over a long
period of time it is found that the average packed was 3 kg with
a standard deviation of 0.1 kg. Assuming the packing is normally
distributed, we can find the fraction (or %) of packages that weigh
more than 3.1 kg.
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">mean</span> <span class="special">=</span> <span class="number">3.</span><span class="special">;</span> <span class="comment">// kg
</span><span class="keyword">double</span> <span class="identifier">standard_deviation</span> <span class="special">=</span> <span class="number">0.1</span><span class="special">;</span> <span class="comment">// kg
</span><span class="identifier">normal</span> <span class="identifier">packs</span><span class="special">(</span><span class="identifier">mean</span><span class="special">,</span> <span class="identifier">standard_deviation</span><span class="special">);</span>
<span class="keyword">double</span> <span class="identifier">max_weight</span> <span class="special">=</span> <span class="number">3.1</span><span class="special">;</span> <span class="comment">// kg
</span><span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Percentage of packs &gt; "</span> <span class="special">&lt;&lt;</span> <span class="identifier">max_weight</span> <span class="special">&lt;&lt;</span> <span class="string">" is "</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">packs</span><span class="special">,</span> <span class="identifier">max_weight</span><span class="special">))</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// P(X &gt; 3.1)
</span>
<span class="keyword">double</span> <span class="identifier">under_weight</span> <span class="special">=</span> <span class="number">2.9</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span><span class="string">"fraction of packs &lt;= "</span> <span class="special">&lt;&lt;</span> <span class="identifier">under_weight</span> <span class="special">&lt;&lt;</span> <span class="string">" with a mean of "</span> <span class="special">&lt;&lt;</span> <span class="identifier">mean</span>
<span class="special">&lt;&lt;</span> <span class="string">" is "</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">packs</span><span class="special">,</span> <span class="identifier">under_weight</span><span class="special">))</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="comment">// fraction of packs &lt;= 2.9 with a mean of 3 is 0.841345
</span><span class="comment">// This is 0.84 - more than the target 0.95
</span><span class="comment">// Want 95% to be over this weight, so what should we set the mean weight to be?
</span><span class="comment">// KK StatCalc says:
</span><span class="keyword">double</span> <span class="identifier">over_mean</span> <span class="special">=</span> <span class="number">3.0664</span><span class="special">;</span>
<span class="identifier">normal</span> <span class="identifier">xpacks</span><span class="special">(</span><span class="identifier">over_mean</span><span class="special">,</span> <span class="identifier">standard_deviation</span><span class="special">);</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"fraction of packs &gt;= "</span> <span class="special">&lt;&lt;</span> <span class="identifier">under_weight</span>
<span class="special">&lt;&lt;</span> <span class="string">" with a mean of "</span> <span class="special">&lt;&lt;</span> <span class="identifier">xpacks</span><span class="special">.</span><span class="identifier">mean</span><span class="special">()</span>
<span class="special">&lt;&lt;</span> <span class="string">" is "</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">xpacks</span><span class="special">,</span> <span class="identifier">under_weight</span><span class="special">))</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="comment">// fraction of packs &gt;= 2.9 with a mean of 3.06449 is 0.950005
</span><span class="keyword">double</span> <span class="identifier">under_fraction</span> <span class="special">=</span> <span class="number">0.05</span><span class="special">;</span> <span class="comment">// so 95% are above the minimum weight mean - sd = 2.9
</span><span class="keyword">double</span> <span class="identifier">low_limit</span> <span class="special">=</span> <span class="identifier">standard_deviation</span><span class="special">;</span>
<span class="keyword">double</span> <span class="identifier">offset</span> <span class="special">=</span> <span class="identifier">mean</span> <span class="special">-</span> <span class="identifier">low_limit</span> <span class="special">-</span> <span class="identifier">quantile</span><span class="special">(</span><span class="identifier">packs</span><span class="special">,</span> <span class="identifier">under_fraction</span><span class="special">);</span>
<span class="keyword">double</span> <span class="identifier">nominal_mean</span> <span class="special">=</span> <span class="identifier">mean</span> <span class="special">+</span> <span class="identifier">offset</span><span class="special">;</span>
<span class="identifier">normal</span> <span class="identifier">nominal_packs</span><span class="special">(</span><span class="identifier">nominal_mean</span><span class="special">,</span> <span class="identifier">standard_deviation</span><span class="special">);</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Setting the packer to "</span> <span class="special">&lt;&lt;</span> <span class="identifier">nominal_mean</span> <span class="special">&lt;&lt;</span> <span class="string">" will mean that "</span>
<span class="special">&lt;&lt;</span> <span class="string">"fraction of packs &gt;= "</span> <span class="special">&lt;&lt;</span> <span class="identifier">under_weight</span>
<span class="special">&lt;&lt;</span> <span class="string">" is "</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">nominal_packs</span><span class="special">,</span> <span class="identifier">under_weight</span><span class="special">))</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span></pre>
<p>
</p>
<p>
</p>
<p>
Setting the packer to 3.06449 will mean that fraction of packs &gt;=
2.9 is 0.95.
</p>
<p>
</p>
<p>
Setting the packer to 3.13263 will mean that fraction of packs &gt;=
2.9 is 0.99, but will more than double the mean loss from 0.0644
to 0.133.
</p>
<p>
</p>
<p>
Alternatively, we could invest in a better (more precise) packer
with a lower standard deviation.
</p>
<p>
</p>
<p>
To estimate how much better (how much smaller standard deviation)
it would have to be, we need to get the 5% quantile to be located
at the under_weight limit, 2.9
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">p</span> <span class="special">=</span> <span class="number">0.05</span><span class="special">;</span> <span class="comment">// wanted p th quantile.
</span><span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Quantile of "</span> <span class="special">&lt;&lt;</span> <span class="identifier">p</span> <span class="special">&lt;&lt;</span> <span class="string">" = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">quantile</span><span class="special">(</span><span class="identifier">packs</span><span class="special">,</span> <span class="identifier">p</span><span class="special">)</span>
<span class="special">&lt;&lt;</span> <span class="string">", mean = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">packs</span><span class="special">.</span><span class="identifier">mean</span><span class="special">()</span> <span class="special">&lt;&lt;</span> <span class="string">", sd = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">packs</span><span class="special">.</span><span class="identifier">standard_deviation</span><span class="special">()</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">//</span></pre>
<p>
</p>
<p>
</p>
<p>
Quantile of 0.05 = 2.83551, mean = 3, sd = 0.1
</p>
<p>
</p>
<p>
With the current packer (mean = 3, sd = 0.1), the 5% quantile is
at 2.8551 kg, a little below our target of 2.9 kg. So we know that
the standard deviation is going to have to be smaller.
</p>
<p>
</p>
<p>
Let's start by guessing that it (now 0.1) needs to be halved, to
a standard deviation of 0.05
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"><span class="identifier">normal</span> <span class="identifier">pack05</span><span class="special">(</span><span class="identifier">mean</span><span class="special">,</span> <span class="number">0.05</span><span class="special">);</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Quantile of "</span> <span class="special">&lt;&lt;</span> <span class="identifier">p</span> <span class="special">&lt;&lt;</span> <span class="string">" = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">quantile</span><span class="special">(</span><span class="identifier">pack05</span><span class="special">,</span> <span class="identifier">p</span><span class="special">)</span>
<span class="special">&lt;&lt;</span> <span class="string">", mean = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">pack05</span><span class="special">.</span><span class="identifier">mean</span><span class="special">()</span> <span class="special">&lt;&lt;</span> <span class="string">", sd = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">pack05</span><span class="special">.</span><span class="identifier">standard_deviation</span><span class="special">()</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span><span class="string">"Fraction of packs &gt;= "</span> <span class="special">&lt;&lt;</span> <span class="identifier">under_weight</span> <span class="special">&lt;&lt;</span> <span class="string">" with a mean of "</span> <span class="special">&lt;&lt;</span> <span class="identifier">mean</span>
<span class="special">&lt;&lt;</span> <span class="string">" and standard deviation of "</span> <span class="special">&lt;&lt;</span> <span class="identifier">pack05</span><span class="special">.</span><span class="identifier">standard_deviation</span><span class="special">()</span>
<span class="special">&lt;&lt;</span> <span class="string">" is "</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">pack05</span><span class="special">,</span> <span class="identifier">under_weight</span><span class="special">))</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="comment">//</span></pre>
<p>
</p>
<p>
</p>
<p>
Fraction of packs &gt;= 2.9 with a mean of 3 and standard deviation
of 0.05 is 0.9772
</p>
<p>
</p>
<p>
So 0.05 was quite a good guess, but we are a little over the 2.9
target, so the standard deviation could be a tiny bit more. So we
could do some more guessing to get closer, say by increasing to 0.06
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"><span class="identifier">normal</span> <span class="identifier">pack06</span><span class="special">(</span><span class="identifier">mean</span><span class="special">,</span> <span class="number">0.06</span><span class="special">);</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Quantile of "</span> <span class="special">&lt;&lt;</span> <span class="identifier">p</span> <span class="special">&lt;&lt;</span> <span class="string">" = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">quantile</span><span class="special">(</span><span class="identifier">pack06</span><span class="special">,</span> <span class="identifier">p</span><span class="special">)</span>
<span class="special">&lt;&lt;</span> <span class="string">", mean = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">pack06</span><span class="special">.</span><span class="identifier">mean</span><span class="special">()</span> <span class="special">&lt;&lt;</span> <span class="string">", sd = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">pack06</span><span class="special">.</span><span class="identifier">standard_deviation</span><span class="special">()</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span><span class="string">"Fraction of packs &gt;= "</span> <span class="special">&lt;&lt;</span> <span class="identifier">under_weight</span> <span class="special">&lt;&lt;</span> <span class="string">" with a mean of "</span> <span class="special">&lt;&lt;</span> <span class="identifier">mean</span>
<span class="special">&lt;&lt;</span> <span class="string">" and standard deviation of "</span> <span class="special">&lt;&lt;</span> <span class="identifier">pack06</span><span class="special">.</span><span class="identifier">standard_deviation</span><span class="special">()</span>
<span class="special">&lt;&lt;</span> <span class="string">" is "</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">pack06</span><span class="special">,</span> <span class="identifier">under_weight</span><span class="special">))</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span></pre>
<p>
</p>
<p>
</p>
<p>
Fraction of packs &gt;= 2.9 with a mean of 3 and standard deviation
of 0.06 is 0.9522
</p>
<p>
</p>
<p>
Now we are getting really close, but to do the job properly, we could
use root finding method, for example the tools provided, and used
elsewhere, in the Math Toolkit, see <a class="link" href="../../../../toolkit/internals1/roots2.html" title="Root Finding Without Derivatives">Root
Finding Without Derivatives</a>.
</p>
<p>
</p>
<p>
But in this normal distribution case, we could be even smarter and
make a direct calculation.
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"><span class="identifier">normal</span> <span class="identifier">s</span><span class="special">;</span> <span class="comment">// For standard normal distribution,
</span><span class="keyword">double</span> <span class="identifier">sd</span> <span class="special">=</span> <span class="number">0.1</span><span class="special">;</span>
<span class="keyword">double</span> <span class="identifier">x</span> <span class="special">=</span> <span class="number">2.9</span><span class="special">;</span> <span class="comment">// Our required limit.
</span><span class="comment">// then probability p = N((x - mean) / sd)
</span><span class="comment">// So if we want to find the standard deviation that would be required to meet this limit,
</span><span class="comment">// so that the p th quantile is located at x,
</span><span class="comment">// in this case the 0.95 (95%) quantile at 2.9 kg pack weight, when the mean is 3 kg.
</span>
<span class="keyword">double</span> <span class="identifier">prob</span> <span class="special">=</span> <span class="identifier">pdf</span><span class="special">(</span><span class="identifier">s</span><span class="special">,</span> <span class="special">(</span><span class="identifier">x</span> <span class="special">-</span> <span class="identifier">mean</span><span class="special">)</span> <span class="special">/</span> <span class="identifier">sd</span><span class="special">);</span>
<span class="keyword">double</span> <span class="identifier">qp</span> <span class="special">=</span> <span class="identifier">quantile</span><span class="special">(</span><span class="identifier">s</span><span class="special">,</span> <span class="number">0.95</span><span class="special">);</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"prob = "</span> <span class="special">&lt;&lt;</span> <span class="identifier">prob</span> <span class="special">&lt;&lt;</span> <span class="string">", quantile(p) "</span> <span class="special">&lt;&lt;</span> <span class="identifier">qp</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// p = 0.241971, quantile(p) 1.64485
</span><span class="comment">// Rearranging, we can directly calculate the required standard deviation:
</span><span class="keyword">double</span> <span class="identifier">sd95</span> <span class="special">=</span> <span class="identifier">abs</span><span class="special">((</span><span class="identifier">x</span> <span class="special">-</span> <span class="identifier">mean</span><span class="special">))</span> <span class="special">/</span> <span class="identifier">qp</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"If we want the "</span><span class="special">&lt;&lt;</span> <span class="identifier">p</span> <span class="special">&lt;&lt;</span> <span class="string">" th quantile to be located at "</span>
<span class="special">&lt;&lt;</span> <span class="identifier">x</span> <span class="special">&lt;&lt;</span> <span class="string">", would need a standard deviation of "</span> <span class="special">&lt;&lt;</span> <span class="identifier">sd95</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">normal</span> <span class="identifier">pack95</span><span class="special">(</span><span class="identifier">mean</span><span class="special">,</span> <span class="identifier">sd95</span><span class="special">);</span> <span class="comment">// Distribution of the 'ideal better' packer.
</span><span class="identifier">cout</span> <span class="special">&lt;&lt;</span><span class="string">"Fraction of packs &gt;= "</span> <span class="special">&lt;&lt;</span> <span class="identifier">under_weight</span> <span class="special">&lt;&lt;</span> <span class="string">" with a mean of "</span> <span class="special">&lt;&lt;</span> <span class="identifier">mean</span>
<span class="special">&lt;&lt;</span> <span class="string">" and standard deviation of "</span> <span class="special">&lt;&lt;</span> <span class="identifier">pack95</span><span class="special">.</span><span class="identifier">standard_deviation</span><span class="special">()</span>
<span class="special">&lt;&lt;</span> <span class="string">" is "</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">pack95</span><span class="special">,</span> <span class="identifier">under_weight</span><span class="special">))</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="comment">// Fraction of packs &gt;= 2.9 with a mean of 3 and standard deviation of 0.0608 is 0.95</span></pre>
<p>
</p>
<p>
</p>
<p>
Notice that these two deceptively simple questions (do we over-fill
or measure better) are actually very common. The weight of beef might
be replaced by a measurement of more or less anything. But the calculations
rely on the accuracy of the standard deviation - something that is
almost always less good than we might wish, especially if based on
a few measurements.
</p>
<p>
</p>
<a name="math_toolkit.dist.stat_tut.weg.normal_example.normal_misc.length_of_bolts"></a><h5>
<a name="id980594"></a>
<a class="link" href="normal_misc.html#math_toolkit.dist.stat_tut.weg.normal_example.normal_misc.length_of_bolts">Length
of bolts</a>
</h5>
<p>
</p>
<p>
A bolt is usable if between 3.9 and 4.1 long. From a large batch
of bolts, a sample of 50 show a mean length of 3.95 with standard
deviation 0.1. Assuming a normal distribution, what proportion is
usable? The true sample mean is unknown, but we can use the sample
mean and standard deviation to find approximate solutions.
</p>
<p>
</p>
<p>
</p>
<pre class="programlisting"> <span class="identifier">normal</span> <span class="identifier">bolts</span><span class="special">(</span><span class="number">3.95</span><span class="special">,</span> <span class="number">0.1</span><span class="special">);</span>
<span class="keyword">double</span> <span class="identifier">top</span> <span class="special">=</span> <span class="number">4.1</span><span class="special">;</span>
<span class="keyword">double</span> <span class="identifier">bottom</span> <span class="special">=</span> <span class="number">3.9</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Fraction long enough [ P(X &lt;= "</span> <span class="special">&lt;&lt;</span> <span class="identifier">top</span> <span class="special">&lt;&lt;</span> <span class="string">") ] is "</span> <span class="special">&lt;&lt;</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">bolts</span><span class="special">,</span> <span class="identifier">top</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Fraction too short [ P(X &lt;= "</span> <span class="special">&lt;&lt;</span> <span class="identifier">bottom</span> <span class="special">&lt;&lt;</span> <span class="string">") ] is "</span> <span class="special">&lt;&lt;</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">bolts</span><span class="special">,</span> <span class="identifier">bottom</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Fraction OK -between "</span> <span class="special">&lt;&lt;</span> <span class="identifier">bottom</span> <span class="special">&lt;&lt;</span> <span class="string">" and "</span> <span class="special">&lt;&lt;</span> <span class="identifier">top</span>
<span class="special">&lt;&lt;</span> <span class="string">"[ P(X &lt;= "</span> <span class="special">&lt;&lt;</span> <span class="identifier">top</span> <span class="special">&lt;&lt;</span> <span class="string">") - P(X&lt;= "</span> <span class="special">&lt;&lt;</span> <span class="identifier">bottom</span> <span class="special">&lt;&lt;</span> <span class="string">" ) ] is "</span>
<span class="special">&lt;&lt;</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">bolts</span><span class="special">,</span> <span class="identifier">top</span><span class="special">)</span> <span class="special">-</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">bolts</span><span class="special">,</span> <span class="identifier">bottom</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"Fraction too long [ P(X &gt; "</span> <span class="special">&lt;&lt;</span> <span class="identifier">top</span> <span class="special">&lt;&lt;</span> <span class="string">") ] is "</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">bolts</span><span class="special">,</span> <span class="identifier">top</span><span class="special">))</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
<span class="identifier">cout</span> <span class="special">&lt;&lt;</span> <span class="string">"95% of bolts are shorter than "</span> <span class="special">&lt;&lt;</span> <span class="identifier">quantile</span><span class="special">(</span><span class="identifier">bolts</span><span class="special">,</span> <span class="number">0.95</span><span class="special">)</span> <span class="special">&lt;&lt;</span> <span class="identifier">endl</span><span class="special">;</span>
</pre>
<p>
</p>
<p>
</p>
</div>
<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
<td align="left"></td>
<td align="right"><div class="copyright-footer">Copyright &#169; 2006 , 2007, 2008, 2009, 2010 John Maddock, Paul A. Bristow,
Hubert Holin, Xiaogang Zhang, Bruno Lalande, Johan R&#229;de, Gautam Sewani and
Thijs van den Berg<p>
Distributed under the Boost Software License, Version 1.0. (See accompanying
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>)
</p>
</div></td>
</tr></table>
<hr>
<div class="spirit-nav">
<a accesskey="p" href="../normal_example.html"><img src="../../../../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../normal_example.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="../nccs_eg.html"><img src="../../../../../../../../../../doc/src/images/next.png" alt="Next"></a>
</div>
</body>
</html>