[boost-doc-zh] r283 commited - 校验，补充和修正
- From: codesite-noreply@xxxxxxxxxx
- To: boost-doc-zh-notify@xxxxxxxxxxxxx
- Date: Wed, 29 Jul 2009 14:44:26 +0000
Revision: 283
Author: alai04
Date: Wed Jul 29 07:28:12 2009
Log: 校验，补充和修正
http://code.google.com/p/boost-doc-zh/source/detail?r=283
Added:
/trunk/doc/html/boost/graph_traits_mpi_graph__id1910405.html
/trunk/doc/html/boost/mpi/allocator_void_id1956361
/trunk/doc/html/boost/mpi/allocator_void_id1956361/rebind.html
/trunk/libs/math/example
/trunk/libs/math/example/Jamfile.v2
/trunk/libs/math/example/binomial_coinflip_example.cpp
/trunk/libs/math/example/binomial_confidence_limits.cpp
/trunk/libs/math/example/binomial_example3.cpp
/trunk/libs/math/example/binomial_example_nag.cpp
/trunk/libs/math/example/binomial_quiz_example.cpp
/trunk/libs/math/example/binomial_sample_sizes.cpp
/trunk/libs/math/example/c_error_policy_example.cpp
/trunk/libs/math/example/chi_square_std_dev_test.cpp
/trunk/libs/math/example/distribution_construction.cpp
/trunk/libs/math/example/error_handling_example.cpp
/trunk/libs/math/example/error_policies_example.cpp
/trunk/libs/math/example/error_policy_example.cpp
/trunk/libs/math/example/f_test.cpp
/trunk/libs/math/example/find_location_example.cpp
/trunk/libs/math/example/find_mean_and_sd_normal.cpp
/trunk/libs/math/example/find_root_example.cpp
/trunk/libs/math/example/find_scale_example.cpp
/trunk/libs/math/example/nc_chi_sq_example.cpp
/trunk/libs/math/example/neg_binom_confidence_limits.cpp
/trunk/libs/math/example/neg_binomial_sample_sizes.cpp
/trunk/libs/math/example/negative_binomial_example1.cpp
/trunk/libs/math/example/negative_binomial_example2.cpp
/trunk/libs/math/example/normal_misc_examples.cpp
/trunk/libs/math/example/policy_eg_1.cpp
/trunk/libs/math/example/policy_eg_10.cpp
/trunk/libs/math/example/policy_eg_2.cpp
/trunk/libs/math/example/policy_eg_3.cpp
/trunk/libs/math/example/policy_eg_4.cpp
/trunk/libs/math/example/policy_eg_5.cpp
/trunk/libs/math/example/policy_eg_6.cpp
/trunk/libs/math/example/policy_eg_7.cpp
/trunk/libs/math/example/policy_eg_8.cpp
/trunk/libs/math/example/policy_eg_9.cpp
/trunk/libs/math/example/policy_ref_snip1.cpp
/trunk/libs/math/example/policy_ref_snip10.cpp
/trunk/libs/math/example/policy_ref_snip11.cpp
/trunk/libs/math/example/policy_ref_snip12.cpp
/trunk/libs/math/example/policy_ref_snip13.cpp
/trunk/libs/math/example/policy_ref_snip2.cpp
/trunk/libs/math/example/policy_ref_snip3.cpp
/trunk/libs/math/example/policy_ref_snip4.cpp
/trunk/libs/math/example/policy_ref_snip5.cpp
/trunk/libs/math/example/policy_ref_snip6.cpp
/trunk/libs/math/example/policy_ref_snip7.cpp
/trunk/libs/math/example/policy_ref_snip8.cpp
/trunk/libs/math/example/policy_ref_snip9.cpp
/trunk/libs/math/example/students_t_example1.cpp
/trunk/libs/math/example/students_t_example2.cpp
/trunk/libs/math/example/students_t_example3.cpp
/trunk/libs/math/example/students_t_single_sample.cpp
/trunk/libs/math/example/students_t_two_samples.cpp
/trunk/libs/math/test/compile_test
/trunk/libs/math/test/compile_test/distribution_concept_check.cpp
/trunk/libs/math/test/std_real_concept_check.cpp
/trunk/libs/math/tools
/trunk/libs/math/tools/lanczos_generator.cpp
/trunk/libs/mpl/example
/trunk/libs/mpl/example/fsm
/trunk/libs/mpl/example/fsm/README.txt
/trunk/libs/mpl/example/fsm/aux_
/trunk/libs/mpl/example/fsm/aux_/STT_impl_gen.hpp
/trunk/libs/mpl/example/fsm/aux_/base_event.hpp
/trunk/libs/mpl/example/fsm/aux_/event.hpp
/trunk/libs/mpl/example/fsm/aux_/state.hpp
/trunk/libs/mpl/example/fsm/aux_/transition.hpp
/trunk/libs/mpl/example/fsm/player.cpp
/trunk/libs/mpl/example/fsm/player1.cpp
/trunk/libs/mpl/example/fsm/player2.cpp
/trunk/libs/mpl/example/fsm/state_machine.hpp
/trunk/libs/mpl/example/inherit_linearly.cpp
/trunk/libs/mpl/example/inherit_multiply.cpp
/trunk/libs/mpl/example/integer.cpp
/trunk/libs/mpl/example/tuple_from_list.cpp
/trunk/libs/python/doc/v2/callbacks.txt
/trunk/libs/python/example/quickstart/Jamroot
/trunk/libs/python/example/quickstart/boost-build.jam
/trunk/libs/python/example/tutorial/Jamroot
/trunk/libs/regex/example
/trunk/libs/regex/example/Jamfile.v2
/trunk/libs/regex/example/grep
/trunk/libs/regex/example/grep/grep.cpp
/trunk/libs/regex/example/snippets
/trunk/libs/regex/example/snippets/captures_example.cpp
/trunk/libs/regex/example/snippets/credit_card_example.cpp
/trunk/libs/regex/example/snippets/icu_example.cpp
/trunk/libs/regex/example/snippets/mfc_example.cpp
/trunk/libs/regex/example/snippets/partial_regex_grep.cpp
/trunk/libs/regex/example/snippets/partial_regex_iterate.cpp
/trunk/libs/regex/example/snippets/partial_regex_match.cpp
/trunk/libs/regex/example/snippets/regex_grep_example_1.cpp
/trunk/libs/regex/example/snippets/regex_grep_example_2.cpp
/trunk/libs/regex/example/snippets/regex_grep_example_3.cpp
/trunk/libs/regex/example/snippets/regex_grep_example_4.cpp
/trunk/libs/regex/example/snippets/regex_iterator_example.cpp
/trunk/libs/regex/example/snippets/regex_match_example.cpp
/trunk/libs/regex/example/snippets/regex_merge_example.cpp
/trunk/libs/regex/example/snippets/regex_replace_example.cpp
/trunk/libs/regex/example/snippets/regex_search_example.cpp
/trunk/libs/regex/example/snippets/regex_split_example_1.cpp
/trunk/libs/regex/example/snippets/regex_split_example_2.cpp
/trunk/libs/regex/example/snippets/regex_token_iterator_eg_1.cpp
/trunk/libs/regex/example/snippets/regex_token_iterator_eg_2.cpp
/trunk/libs/regex/example/timer
/trunk/libs/regex/example/timer/bc55.mak
/trunk/libs/regex/example/timer/bcb4.mak
/trunk/libs/regex/example/timer/bcb5.mak
/trunk/libs/regex/example/timer/gcc.mak
/trunk/libs/regex/example/timer/input_script.txt
/trunk/libs/regex/example/timer/regex_timer.cpp
/trunk/libs/regex/example/timer/vc6-stlport.mak
/trunk/libs/regex/example/timer/vc6.mak
/trunk/libs/wave/doc/class_ref_ctxpolicy_depr.html
Deleted:
/trunk/doc/html/BOOST_PP_LOCAL_LIMITS_id2868241.html
/trunk/doc/html/BOOST_PP_LOCAL_LIMITS_id2872702.html
/trunk/doc/html/BOOST_PP_LOCAL_LIMITS_id2876989.html
/trunk/doc/html/BOOST_PP_LOCAL_LIMITS_id2880111.html
/trunk/doc/html/BOOST_PP_LOCAL_LIMITS_id2884330.html
/trunk/doc/html/BOOST_PP_LOCAL_LIMITS_id2888376.html
/trunk/doc/html/BOOST_PP_LOCAL_LIMITS_id2891852.html
/trunk/doc/html/BOOST_PP_LOCAL_LIMITS_id2899680.html
/trunk/doc/html/BOOST_PP_LOCAL_MACRO_id2868232.html
/trunk/doc/html/BOOST_PP_LOCAL_MACRO_id2872692.html
/trunk/doc/html/BOOST_PP_LOCAL_MACRO_id2876980.html
/trunk/doc/html/BOOST_PP_LOCAL_MACRO_id2880102.html
/trunk/doc/html/BOOST_PP_LOCAL_MACRO_id2884320.html
/trunk/doc/html/BOOST_PP_LOCAL_MACRO_id2888366.html
/trunk/doc/html/BOOST_PP_LOCAL_MACRO_id2891842.html
/trunk/doc/html/BOOST_PP_LOCAL_MACRO_id2899671.html
/trunk/doc/html/boost/graph_traits_mpi_graph__id2665953.html
/trunk/doc/html/boost/mpi/allocator_void_id2813114
Modified:
/trunk/doc/html/boost_tr1/header_list.html
/trunk/doc/html/boost_tr1/subject_list.html
/trunk/doc/html/date_time/doxy.html
/trunk/doc/html/libraries.html
/trunk/doc/html/program_options/design.html
/trunk/doc/html/program_options/reference.html
/trunk/doc/html/program_options/tutorial.html
/trunk/doc/html/proto/reference.html
/trunk/doc/html/thread/changes.html
/trunk/doc/html/thread/synchronization.html
/trunk/doc/html/thread/thread_local_storage.html
/trunk/doc/html/thread/thread_management.html
/trunk/doc/html/thread/time.html
/trunk/libs/exception/doc/error_info.html
/trunk/libs/exception/doc/tutorial_transporting_data.html
/trunk/libs/graph/test/adj_list_cc.cpp
/trunk/libs/graph/test/adj_list_edge_list_set.cpp
/trunk/libs/graph/test/adj_list_test.cpp
/trunk/libs/graph/test/adj_matrix_cc.cpp
/trunk/libs/graph/test/adjacency_matrix_test.cpp
/trunk/libs/graph/test/all_planar_input_files_test.cpp
/trunk/libs/graph/test/astar_search_test.cpp
/trunk/libs/graph/test/basic_planarity_test.cpp
/trunk/libs/graph/test/bellman-test.cpp
/trunk/libs/graph/test/betweenness_centrality_test.cpp
/trunk/libs/graph/test/bfs.cpp
/trunk/libs/graph/test/bfs_cc.cpp
/trunk/libs/graph/test/biconnected_components_test.cpp
/trunk/libs/graph/test/bidir_remove_edge.cpp
/trunk/libs/graph/test/bidir_vec_remove_edge.cpp
/trunk/libs/graph/test/bundled_properties.cpp
/trunk/libs/graph/test/copy.cpp
/trunk/libs/graph/test/csr_graph_test.cpp
/trunk/libs/graph/test/cuthill_mckee_ordering.cpp
/trunk/libs/graph/test/cycle_ratio_tests.cpp
/trunk/libs/graph/test/dag_longest_paths.cpp
/trunk/libs/graph/test/dfs.cpp
/trunk/libs/graph/test/dfs_cc.cpp
/trunk/libs/graph/test/dijkstra_cc.cpp
/trunk/libs/graph/test/dijkstra_heap_performance.cpp
/trunk/libs/graph/test/dominator_tree_test.cpp
/trunk/libs/graph/test/edge_list_cc.cpp
/trunk/libs/graph/test/filter_graph_vp_test.cpp
/trunk/libs/graph/test/filtered_graph_cc.cpp
/trunk/libs/graph/test/floyd_warshall_test.cpp
/trunk/libs/graph/test/graph.cpp
/trunk/libs/graph/test/graph_concepts.cpp
/trunk/libs/graph/test/graph_type.hpp
/trunk/libs/graph/test/graphml_test.cpp
/trunk/libs/graph/test/graphml_test.xml
/trunk/libs/graph/test/graphviz_test.cpp
/trunk/libs/graph/test/gursoy_atun_layout_test.cpp
/trunk/libs/graph/test/isomorphism.cpp
/trunk/libs/graph/test/johnson-test.cpp
/trunk/libs/graph/test/king_ordering.cpp
/trunk/libs/graph/test/kolmogorov_max_flow_test.cpp
/trunk/libs/graph/test/layout_test.cpp
/trunk/libs/graph/test/leda_graph_cc.cpp
/trunk/libs/graph/test/lvalue_pmap.cpp
/trunk/libs/graph/test/make_bicon_planar_test.cpp
/trunk/libs/graph/test/make_connected_test.cpp
/trunk/libs/graph/test/make_maximal_planar_test.cpp
/trunk/libs/graph/test/matching_test.cpp
/trunk/libs/graph/test/max_flow_test.cpp
/trunk/libs/graph/test/parallel_edges_loops_test.cpp
/trunk/libs/graph/test/property_iter.cpp
/trunk/libs/graph/test/random_matching_test.cpp
/trunk/libs/graph/test/relaxed_heap_test.cpp
/trunk/libs/graph/test/reverse_graph_cc.cpp
/trunk/libs/graph/test/sequential_vertex_coloring.cpp
/trunk/libs/graph/test/serialize.cpp
/trunk/libs/graph/test/stanford_graph_cc.cpp
/trunk/libs/graph/test/subgraph.cpp
/trunk/libs/graph/test/transitive_closure_test.cpp
/trunk/libs/graph/test/transitive_closure_test2.cpp
/trunk/libs/graph/test/vector_graph_cc.cpp
/trunk/libs/math/doc/sf_and_dist/html/math_toolkit/main_overview/tr1.html
/trunk/libs/type_traits/doc/html/boost_typetraits/reference/add_cv.html
/trunk/libs/type_traits/doc/html/boost_typetraits/reference/add_pointer.html
/trunk/libs/type_traits/doc/html/boost_typetraits/reference/add_reference.html
/trunk/libs/type_traits/doc/html/boost_typetraits/reference/add_volatile.html
=======================================
--- /dev/null
+++ /trunk/doc/html/boost/graph_traits_mpi_graph__id1910405.html Wed Jul 29
07:28:12 2009
@@ -0,0 +1,78 @@
+<html>
+<head>
+<meta http-equiv="Content-Type" content="text/html; charset=ISO-8859-1">
+<title>Struct graph_traits<mpi::graph_communicator></title>
+<link rel="stylesheet" href="../boostbook.css" type="text/css">
+<meta name="generator" content="DocBook XSL Stylesheets V1.74.3">
+<link rel="home" href="../index.html" title="The Boost C++ Libraries
BoostBook Documentation Subset">
+<link rel="up"
href="../mpi/reference.html#header.boost.mpi.graph_communicator_hpp"
title="Header <boost/mpi/graph_communicator.hpp>">
+<link rel="prev" href="mpi/get_id1911065.html" title="Function get">
+<link rel="next" href="mpi/group.html" title="Class group">
+</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="mpi/get_id1911065.html"><img
src="../../../doc/html/images/prev.png" alt="Prev"></a><a accesskey="u"
href="../mpi/reference.html#header.boost.mpi.graph_communicator_hpp"><img
src="../../../doc/html/images/up.png" alt="Up"></a><a accesskey="h"
href="../index.html"><img src="../../../doc/html/images/home.png"
alt="Home"></a><a accesskey="n" href="mpi/group.html"><img
src="../../../doc/html/images/next.png" alt="Next"></a>
+</div>
+<div class="refentry" lang="en">
+<a name="boost.graph_traits_mpi_graph__id1910405"></a><div
class="titlepage"></div>
+<div class="refnamediv">
+<h2><span class="refentrytitle">Struct
graph_traits<mpi::graph_communicator></span></h2>
+<p>boost::graph_traits<mpi::graph_communicator> — Traits
structure that allows a communicator with graph topology to be view as a
graph by the Boost Graph Library. </p>
+</div>
+<h2 xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision";
class="refsynopsisdiv-title">Synopsis</h2>
+<div xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision";
class="refsynopsisdiv"><pre class="synopsis"><span class="emphasis"><em>//
In header: <<a class="link"
href="../mpi/reference.html#header.boost.mpi.graph_communicator_hpp"
title="Header
<boost/mpi/graph_communicator.hpp>">boost/mpi/graph_communicator.hpp</a>>
+
+</em></span>
+<span class="bold"><strong>struct</strong></span> <a class="link"
href="graph_traits_mpi_graph__id1910405.html" title="Struct
graph_traits<mpi::graph_communicator>">graph_traits</a><mpi::graph_communicator>
{
+ <span class="emphasis"><em>// types</em></span>
+ <span class="bold"><strong>typedef</strong></span> <span
class="bold"><strong>int</strong></span> <a
name="boost.graph_traits_mpi_graph__id1910405.vertex_descriptor"></a>vertex_descriptor;
+ <span class="bold"><strong>typedef</strong></span> std::pair< <span
class="bold"><strong>int</strong></span>, <span
class="bold"><strong>int</strong></span> > <a
name="boost.graph_traits_mpi_graph__id1910405.edge_descriptor"></a>edge_descriptor;
+ <span class="bold"><strong>typedef</strong></span>
directed_tag <a
name="boost.graph_traits_mpi_graph__id1910405.directed_category"></a>directed_category;
+ <span class="bold"><strong>typedef</strong></span>
disallow_parallel_edge_tag <a
name="boost.graph_traits_mpi_graph__id1910405.edge_parallel_category"></a>edge_parallel_category;
+ <span class="bold"><strong>typedef</strong></span> <span
class="emphasis"><em>unspecified</em></span> <a
name="boost.graph_traits_mpi_graph__id1910405.out_edge_iterator"></a>out_edge_iterator;
+ <span class="bold"><strong>typedef</strong></span> <span
class="bold"><strong>int</strong></span> <a
name="boost.graph_traits_mpi_graph__id1910405.degree_size_type"></a>degree_size_type;
+ <span class="bold"><strong>typedef</strong></span> <span
class="emphasis"><em>unspecified</em></span> <a
name="boost.graph_traits_mpi_graph__id1910405.adjacency_iterator"></a>adjacency_iterator;
+ <span class="bold"><strong>typedef</strong></span> counting_iterator<
<span class="bold"><strong>int</strong></span> > <a
name="boost.graph_traits_mpi_graph__id1910405.vertex_iterator"></a>vertex_iterator;
+ <span class="bold"><strong>typedef</strong></span> <span
class="bold"><strong>int</strong></span> <a
name="boost.graph_traits_mpi_graph__id1910405.vertices_size_type"></a>vertices_size_type;
+ <span class="bold"><strong>typedef</strong></span> <span
class="emphasis"><em>unspecified</em></span> <a
name="boost.graph_traits_mpi_graph__id1910405.edge_iterator"></a>edge_iterator;
+ <span class="bold"><strong>typedef</strong></span> <span
class="bold"><strong>int</strong></span> <a
name="boost.graph_traits_mpi_graph__id1910405.edges_size_type"></a>edges_size_type;
+
+ <span class="emphasis"><em>// <a class="link"
href="graph_traits_mpi_graph__id1910405.html#id1910498-bb">public static
functions</a></em></span>
+ <span class="type"><span class="bold"><strong>static</strong></span>
vertex_descriptor</span> <a class="link"
href="graph_traits_mpi_graph__id1910405.html#id1910502-bb">null_vertex</a>() ;
+};</pre></div>
+<div class="refsect1" lang="en">
+<a name="id2860204"></a><h2>Description</h2>
+<p>The specialization of <code class="computeroutput">graph_traits</code>
for an MPI communicator allows a communicator with graph topology to be
viewed as a graph. An MPI communicator with graph topology meets the
requirements of the Graph, Incidence Graph, Adjacency Graph, Vertex List
Graph, and Edge List Graph concepts from the Boost Graph Library. </p>
+<div class="refsect2" lang="en">
+<a name="id2860221"></a><h3>
+<a name="id1910498-bb"></a><code
class="computeroutput">graph_traits</code> public static functions</h3>
+<div class="orderedlist"><ol type="1"><li>
+<pre class="literallayout"><span class="type"><span
class="bold"><strong>static</strong></span> vertex_descriptor</span> <a
name="id1910502-bb"></a>null_vertex() ;</pre>Returns a vertex descriptor
that can never refer to any valid vertex. </li></ol></div>
+</div>
+</div>
+</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 © 2005-2007
Douglas Gregor,
+ Matthias Troyer, Trustees of Indiana University<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="mpi/get_id1911065.html"><img
src="../../../doc/html/images/prev.png" alt="Prev"></a><a accesskey="u"
href="../mpi/reference.html#header.boost.mpi.graph_communicator_hpp"><img
src="../../../doc/html/images/up.png" alt="Up"></a><a accesskey="h"
href="../index.html"><img src="../../../doc/html/images/home.png"
alt="Home"></a><a accesskey="n" href="mpi/group.html"><img
src="../../../doc/html/images/next.png" alt="Next"></a>
+</div>
+</body>
+</html>
=======================================
--- /dev/null
+++ /trunk/doc/html/boost/mpi/allocator_void_id1956361/rebind.html Wed Jul
29 07:28:12 2009
@@ -0,0 +1,56 @@
+<html>
+<head>
+<meta http-equiv="Content-Type" content="text/html; charset=ISO-8859-1">
+<title>Struct template rebind</title>
+<link rel="stylesheet" href="../../../boostbook.css" type="text/css">
+<meta name="generator" content="DocBook XSL Stylesheets V1.74.3">
+<link rel="home" href="../../../index.html" title="The Boost C++ Libraries
BoostBook Documentation Subset">
+<link rel="up" href="../allocator_void_id1956361.html#id2825395"
title="Description">
+<link rel="prev" href="../../../mpi/reference.html" title="Reference">
+<link rel="next" href="../allocator.html" title="Class template allocator">
+</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="../../../mpi/reference.html"><img
src="../../../../../doc/html/images/prev.png" alt="Prev"></a><a
accesskey="u" href="../allocator_void_id1956361.html#id2825395"><img
src="../../../../../doc/html/images/up.png" alt="Up"></a><a accesskey="h"
href="../../../index.html"><img
src="../../../../../doc/html/images/home.png" alt="Home"></a><a
accesskey="n" href="../allocator.html"><img
src="../../../../../doc/html/images/next.png" alt="Next"></a>
+</div>
+<div class="refentry" lang="en">
+<a name="boost.mpi.allocator_void_id1956361.rebind"></a><div
class="titlepage"></div>
+<div class="refnamediv">
+<h2><span class="refentrytitle">Struct template rebind</span></h2>
+<p>boost::mpi::allocator<void>::rebind</p>
+</div>
+<h2 xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision";
class="refsynopsisdiv-title">Synopsis</h2>
+<div xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision";
class="refsynopsisdiv"><pre class="synopsis"><span class="emphasis"><em>//
In header: <<a class="link"
href="../../../mpi/reference.html#header.boost.mpi.allocator_hpp"
title="Header
<boost/mpi/allocator.hpp>">boost/mpi/allocator.hpp</a>>
+
+</em></span>
+<span class="bold"><strong>template</strong></span><<span
class="bold"><strong>typename</strong></span> U>
+<span class="bold"><strong>struct</strong></span> <a class="link"
href="rebind.html" title="Struct template rebind">rebind</a> {
+ <span class="emphasis"><em>// types</em></span>
+ <span class="bold"><strong>typedef</strong></span> <a class="link"
href="../allocator.html" title="Class template allocator">allocator</a><
U > <a name="boost.mpi.allocator_void_id1956361.rebind.other"></a>other;
+};</pre></div>
+</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 © 2005-2007
Douglas Gregor,
+ Matthias Troyer, Trustees of Indiana University<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="../../../mpi/reference.html"><img
src="../../../../../doc/html/images/prev.png" alt="Prev"></a><a
accesskey="u" href="../allocator_void_id1956361.html#id2825395"><img
src="../../../../../doc/html/images/up.png" alt="Up"></a><a accesskey="h"
href="../../../index.html"><img
src="../../../../../doc/html/images/home.png" alt="Home"></a><a
accesskey="n" href="../allocator.html"><img
src="../../../../../doc/html/images/next.png" alt="Next"></a>
+</div>
+</body>
+</html>
=======================================
--- /dev/null
+++ /trunk/libs/math/example/Jamfile.v2 Wed Jul 29 07:28:12 2009
@@ -0,0 +1,79 @@
+# \math_toolkit\libs\math\example\jamfile.v2
+# Runs statistics examples
+# Copyright 2007 John Maddock and Paul A. Bristow.
+# Distributed under the Boost Software License, Version 1.0.
+# (See accompanying file LICENSE_1_0.txt or copy at
http://www.boost.org/LICENSE_1_0.txt)
+
+# bring in the rules for testing
+import testing ;
+
+project
+ : requirements
+ <toolset>gcc:<cxxflags>-Wno-missing-braces
+ <toolset>darwin:<cxxflags>-Wno-missing-braces
+ <toolset>acc:<cxxflags>+W2068,2461,2236,4070
+ <toolset>intel:<cxxflags>-Qwd264,239
+ <toolset>msvc:<warnings>all
+ <toolset>msvc:<asynch-exceptions>on
+ <toolset>msvc:<cxxflags>/wd4996
+ <toolset>msvc:<cxxflags>/wd4512
+ <toolset>msvc:<cxxflags>/wd4610
+ <toolset>msvc:<cxxflags>/wd4510
+ <toolset>msvc:<cxxflags>/wd4127
+ <toolset>msvc:<cxxflags>/wd4701
+ <toolset>msvc:<cxxflags>/wd4127
+ <toolset>msvc:<cxxflags>/wd4305
+ <include>../../..
+ ;
+
+run binomial_confidence_limits.cpp ;
+run binomial_example3.cpp ;
+run binomial_sample_sizes.cpp ;
+run binomial_example_nag.cpp ;
+run c_error_policy_example.cpp ;
+run chi_square_std_dev_test.cpp ;
+run distribution_construction.cpp ;
+run error_handling_example.cpp ;
+run error_policies_example.cpp ;
+run error_policy_example.cpp ;
+run f_test.cpp ;
+run neg_binom_confidence_limits.cpp ;
+run neg_binomial_sample_sizes.cpp ;
+# run negative_binomial_construction_examples.cpp ; delete for now
+run negative_binomial_example1.cpp ;
+run negative_binomial_example2.cpp ;
+# run negative_binomial_example3.cpp ;
+run policy_eg_1.cpp ;
+run policy_eg_2.cpp ;
+run policy_eg_3.cpp ;
+run policy_eg_4.cpp ;
+run policy_eg_5.cpp ;
+run policy_eg_6.cpp ;
+run policy_eg_7.cpp ;
+run policy_eg_8.cpp ;
+run policy_eg_9.cpp ;
+run policy_eg_10.cpp ;
+run policy_ref_snip1.cpp ;
+run policy_ref_snip10.cpp ;
+run policy_ref_snip11.cpp ;
+run policy_ref_snip12.cpp ;
+run policy_ref_snip13.cpp ;
+run policy_ref_snip2.cpp ;
+run policy_ref_snip3.cpp ;
+run policy_ref_snip4.cpp ;
+run policy_ref_snip5.cpp ;
+run policy_ref_snip6.cpp ;
+run policy_ref_snip7.cpp ;
+run policy_ref_snip8.cpp ;
+run policy_ref_snip9.cpp ;
+#run statistics_functions_example1.cpp ;
+run students_t_example1.cpp ;
+run students_t_example2.cpp ;
+run students_t_example3.cpp ;
+run students_t_single_sample.cpp ;
+run students_t_two_samples.cpp ;
+
+
+
+
+
=======================================
--- /dev/null
+++ /trunk/libs/math/example/binomial_coinflip_example.cpp Wed Jul 29
07:28:12 2009
@@ -0,0 +1,237 @@
+// Copyright Paul A. 2007
+// Copyright John Maddock 2006
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Simple example of computing probabilities and quantiles for
+// a Bernoulli random variable representing the flipping of a coin.
+
+// http://mathworld.wolfram.com/CoinTossing.html
+// http://en.wikipedia.org/wiki/Bernoulli_trial
+// Weisstein, Eric W. "Dice." From MathWorld--A Wolfram Web Resource.
+// http://mathworld.wolfram.com/Dice.html
+// http://en.wikipedia.org/wiki/Bernoulli_distribution
+// http://mathworld.wolfram.com/BernoulliDistribution.html
+//
+// An idealized coin consists of a circular disk of zero thickness which,
+// when thrown in the air and allowed to fall, will rest with either side
face up
+// ("heads" H or "tails" T) with equal probability. A coin is therefore a
two-sided die.
+// Despite slight differences between the sides and nonzero thickness of
actual coins,
+// the distribution of their tosses makes a good approximation to a p==1/2
Bernoulli distribution.
+
+//[binomial_coinflip_example1
+
+/*`An example of a [@http://en.wikipedia.org/wiki/Bernoulli_process
Bernoulli process]
+is coin flipping.
+A variable in such a sequence may be called a Bernoulli variable.
+
+This example shows using the Binomial distribution to predict the
probability
+of heads and tails when throwing a coin.
+
+The number of correct answers (say heads),
+X, is distributed as a binomial random variable
+with binomial distribution parameters number of trials (flips) n = 10 and
probability (success_fraction) of getting a head p = 0.5 (a 'fair' coin).
+
+(Our coin is assumed fair, but we could easily change the success_fraction
parameter p
+from 0.5 to some other value to simulate an unfair coin,
+say 0.6 for one with chewing gum on the tail,
+so it is more likely to fall tails down and heads up).
+
+First we need some includes and using statements to be able to use the
binomial distribution, some std input and output, and get started:
+*/
+
+#include <boost/math/distributions/binomial.hpp>
+ using boost::math::binomial;
+
+#include <iostream>
+ using std::cout; using std::endl; using std::left;
+#include <iomanip>
+ using std::setw;
+
+int main()
+{
+ cout << "Using Binomial distribution to predict how many heads and
tails." << endl;
+ try
+ {
+/*`
+See note [link coinflip_eg_catch with the catch block]
+about why a try and catch block is always a good idea.
+
+First, construct a binomial distribution with parameters success_fraction
+1/2, and how many flips.
+*/
+ const double success_fraction = 0.5; // = 50% = 1/2 for a 'fair' coin.
+ int flips = 10;
+ binomial flip(flips, success_fraction);
+
+ cout.precision(4);
+/*`
+ Then some examples of using Binomial moments (and echoing the parameters).
+*/
+ cout << "From " << flips << " one can expect to get on average "
+ << mean(flip) << " heads (or tails)." << endl;
+ cout << "Mode is " << mode(flip) << endl;
+ cout << "Standard deviation is " << standard_deviation(flip) << endl;
+ cout << "So about 2/3 will lie within 1 standard deviation and get
between "
+ << ceil(mean(flip) - standard_deviation(flip)) << " and "
+ << floor(mean(flip) + standard_deviation(flip)) << " correct." <<
endl;
+ cout << "Skewness is " << skewness(flip) << endl;
+ // Skewness of binomial distributions is only zero (symmetrical)
+ // if success_fraction is exactly one half,
+ // for example, when flipping 'fair' coins.
+ cout << "Skewness if success_fraction is " << flip.success_fraction()
+ << " is " << skewness(flip) << endl << endl; // Expect zero for
a 'fair' coin.
+/*`
+Now we show a variety of predictions on the probability of heads:
+*/
+ cout << "For " << flip.trials() << " coin flips: " << endl;
+ cout << "Probability of getting no heads is " << pdf(flip, 0) << endl;
+ cout << "Probability of getting at least one head is " << 1. -
pdf(flip, 0) << endl;
+/*`
+When we want to calculate the probability for a range or values we can sum
the PDF's:
+*/
+ cout << "Probability of getting 0 or 1 heads is "
+ << pdf(flip, 0) + pdf(flip, 1) << endl; // sum of exactly ==
probabilities
+/*`
+Or we can use the cdf.
+*/
+ cout << "Probability of getting 0 or 1 (<= 1) heads is " << cdf(flip,
1) << endl;
+ cout << "Probability of getting 9 or 10 heads is " << pdf(flip, 9) +
pdf(flip, 10) << endl;
+/*`
+Note that using
+*/
+ cout << "Probability of getting 9 or 10 heads is " << 1. - cdf(flip,
8) << endl;
+/*`
+is less accurate than using the complement
+*/
+ cout << "Probability of getting 9 or 10 heads is " <<
cdf(complement(flip, 8)) << endl;
+/*`
+Since the subtraction may involve
+[@http://docs.sun.com/source/806-3568/ncg_goldberg.html cancellation
error],
+where as `cdf(complement(flip, 8))`
+does not use such a subtraction internally, and so does not exhibit the
problem.
+
+To get the probability for a range of heads, we can either add the pdfs
for each number of heads
+*/
+ cout << "Probability of between 4 and 6 heads (4 or 5 or 6) is "
+ // P(X == 4) + P(X == 5) + P(X == 6)
+ << pdf(flip, 4) + pdf(flip, 5) + pdf(flip, 6) << endl;
+/*`
+But this is probably less efficient than using the cdf
+*/
+ cout << "Probability of between 4 and 6 heads (4 or 5 or 6) is "
+ // P(X <= 6) - P(X <= 3) == P(X < 4)
+ << cdf(flip, 6) - cdf(flip, 3) << endl;
+/*`
+Certainly for a bigger range like, 3 to 7
+*/
+ cout << "Probability of between 3 and 7 heads (3, 4, 5, 6 or 7) is "
+ // P(X <= 7) - P(X <= 2) == P(X < 3)
+ << cdf(flip, 7) - cdf(flip, 2) << endl;
+ cout << endl;
+
+/*`
+Finally, print two tables of probability for the /exactly/ and /at least/
a number of heads.
+*/
+ // Print a table of probability for the exactly a number of heads.
+ cout << "Probability of getting exactly (==) heads" << endl;
+ for (int successes = 0; successes <= flips; successes++)
+ { // Say success means getting a head (or equally success means
getting a tail).
+ double probability = pdf(flip, successes);
+ cout << left << setw(2) << successes << " " << setw(10)
+ << probability << " or 1 in " << 1. / probability
+ << ", or " << probability * 100. << "%" << endl;
+ } // for i
+ cout << endl;
+
+ // Tabulate the probability of getting between zero heads and 0 upto
10 heads.
+ cout << "Probability of getting upto (<=) heads" << endl;
+ for (int successes = 0; successes <= flips; successes++)
+ { // Say success means getting a head
+ // (equally success could mean getting a tail).
+ double probability = cdf(flip, successes); // P(X <= heads)
+ cout << setw(2) << successes << " " << setw(10) << left
+ << probability << " or 1 in " << 1. / probability << ", or "
+ << probability * 100. << "%"<< endl;
+ } // for i
+/*`
+The last (0 to 10 heads) must, of course, be 100% probability.
+*/
+ }
+ catch(const std::exception& e)
+ {
+ //
+ /*`
+ [#coinflip_eg_catch]
+ It is always essential to include try & catch blocks because
+ default policies are to throw exceptions on arguments that
+ are out of domain or cause errors like numeric-overflow.
+
+ Lacking try & catch blocks, the program will abort, whereas the
+ message below from the thrown exception will give some helpful
+ clues as to the cause of the problem.
+ */
+ std::cout <<
+ "\n""Message from thrown exception was:\n " << e.what() <<
std::endl;
+ }
+//] [binomial_coinflip_example1]
+ return 0;
+} // int main()
+
+// Output:
+
+//[binomial_coinflip_example_output
+/*`
+
+[pre
+Using Binomial distribution to predict how many heads and tails.
+From 10 one can expect to get on average 5 heads (or tails).
+Mode is 5
+Standard deviation is 1.581
+So about 2/3 will lie within 1 standard deviation and get between 4 and 6
correct.
+Skewness is 0
+Skewness if success_fraction is 0.5 is 0
+
+For 10 coin flips:
+Probability of getting no heads is 0.0009766
+Probability of getting at least one head is 0.999
+Probability of getting 0 or 1 heads is 0.01074
+Probability of getting 0 or 1 (<= 1) heads is 0.01074
+Probability of getting 9 or 10 heads is 0.01074
+Probability of getting 9 or 10 heads is 0.01074
+Probability of getting 9 or 10 heads is 0.01074
+Probability of between 4 and 6 heads (4 or 5 or 6) is 0.6562
+Probability of between 4 and 6 heads (4 or 5 or 6) is 0.6563
+Probability of between 3 and 7 heads (3, 4, 5, 6 or 7) is 0.8906
+
+Probability of getting exactly (==) heads
+0 0.0009766 or 1 in 1024, or 0.09766%
+1 0.009766 or 1 in 102.4, or 0.9766%
+2 0.04395 or 1 in 22.76, or 4.395%
+3 0.1172 or 1 in 8.533, or 11.72%
+4 0.2051 or 1 in 4.876, or 20.51%
+5 0.2461 or 1 in 4.063, or 24.61%
+6 0.2051 or 1 in 4.876, or 20.51%
+7 0.1172 or 1 in 8.533, or 11.72%
+8 0.04395 or 1 in 22.76, or 4.395%
+9 0.009766 or 1 in 102.4, or 0.9766%
+10 0.0009766 or 1 in 1024, or 0.09766%
+
+Probability of getting upto (<=) heads
+0 0.0009766 or 1 in 1024, or 0.09766%
+1 0.01074 or 1 in 93.09, or 1.074%
+2 0.05469 or 1 in 18.29, or 5.469%
+3 0.1719 or 1 in 5.818, or 17.19%
+4 0.377 or 1 in 2.653, or 37.7%
+5 0.623 or 1 in 1.605, or 62.3%
+6 0.8281 or 1 in 1.208, or 82.81%
+7 0.9453 or 1 in 1.058, or 94.53%
+8 0.9893 or 1 in 1.011, or 98.93%
+9 0.999 or 1 in 1.001, or 99.9%
+10 1 or 1 in 1, or 100%
+]
+*/
+//][/binomial_coinflip_example_output]
=======================================
--- /dev/null
+++ /trunk/libs/math/example/binomial_confidence_limits.cpp Wed Jul 29
07:28:12 2009
@@ -0,0 +1,159 @@
+// Copyright John Maddock 2006
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifdef _MSC_VER
+# pragma warning(disable: 4512) // assignment operator could not be
generated.
+# pragma warning(disable: 4510) // default constructor could not be
generated.
+# pragma warning(disable: 4610) // can never be instantiated - user
defined constructor required.
+#endif
+
+#include <iostream>
+#include <iomanip>
+#include <boost/math/distributions/binomial.hpp>
+
+void confidence_limits_on_frequency(unsigned trials, unsigned successes)
+{
+ //
+ // trials = Total number of trials.
+ // successes = Total number of observed successes.
+ //
+ // Calculate confidence limits for an observed
+ // frequency of occurrence that follows a binomial
+ // distribution.
+ //
+ using namespace std;
+ using namespace boost::math;
+
+ // Print out general info:
+ cout <<
+ "___________________________________________\n"
+ "2-Sided Confidence Limits For Success Ratio\n"
+ "___________________________________________\n\n";
+ cout << setprecision(7);
+ cout << setw(40) << left << "Number of Observations" << "= " << trials
<< "\n";
+ cout << setw(40) << left << "Number of successes" << "= " << successes
<< "\n";
+ cout << setw(40) << left << "Sample frequency of occurrence" << "= "
<< double(successes) / trials << "\n";
+ //
+ // Define a table of significance levels:
+ //
+ double alpha[] = { 0.5, 0.25, 0.1, 0.05, 0.01, 0.001, 0.0001, 0.00001 };
+ //
+ // Print table header:
+ //
+ cout << "\n\n"
+
"_______________________________________________________________________\n"
+ "Confidence Lower CP Upper CP Lower JP
Upper JP\n"
+ " Value (%) Limit Limit Limit
Limit\n"
+
"_______________________________________________________________________\n";
+ //
+ // Now print out the data for the table rows.
+ //
+ for(unsigned i = 0; i < sizeof(alpha)/sizeof(alpha[0]); ++i)
+ {
+ // Confidence value:
+ cout << fixed << setprecision(3) << setw(10) << right << 100 *
(1-alpha[i]);
+ // Calculate Clopper Pearson bounds:
+ double l = binomial_distribution<>::find_lower_bound_on_p(trials,
successes, alpha[i]/2);
+ double u = binomial_distribution<>::find_upper_bound_on_p(trials,
successes, alpha[i]/2);
+ // Print Clopper Pearson Limits:
+ cout << fixed << setprecision(5) << setw(15) << right << l;
+ cout << fixed << setprecision(5) << setw(15) << right << u;
+ // Calculate Jeffreys Prior Bounds:
+ l = binomial_distribution<>::find_lower_bound_on_p(trials,
successes, alpha[i]/2, binomial_distribution<>::jeffreys_prior_interval);
+ u = binomial_distribution<>::find_upper_bound_on_p(trials,
successes, alpha[i]/2, binomial_distribution<>::jeffreys_prior_interval);
+ // Print Jeffreys Prior Limits:
+ cout << fixed << setprecision(5) << setw(15) << right << l;
+ cout << fixed << setprecision(5) << setw(15) << right << u <<
std::endl;
+ }
+ cout << endl;
+} // void confidence_limits_on_frequency()
+
+int main()
+{
+ confidence_limits_on_frequency(20, 4);
+ confidence_limits_on_frequency(200, 40);
+ confidence_limits_on_frequency(2000, 400);
+
+ return 0;
+} // int main()
+
+/*
+
+------ Build started: Project: binomial_confidence_limits, Configuration:
Debug Win32 ------
+Compiling...
+binomial_confidence_limits.cpp
+Linking...
+Autorun
"i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\binomial_confidence_limits.exe"
+___________________________________________
+2-Sided Confidence Limits For Success Ratio
+___________________________________________
+
+Number of Observations = 20
+Number of successes = 4
+Sample frequency of occurrence = 0.2
+
+
+_______________________________________________________________________
+Confidence Lower CP Upper CP Lower JP Upper JP
+ Value (%) Limit Limit Limit Limit
+_______________________________________________________________________
+ 50.000 0.12840 0.29588 0.14974 0.26916
+ 75.000 0.09775 0.34633 0.11653 0.31861
+ 90.000 0.07135 0.40103 0.08734 0.37274
+ 95.000 0.05733 0.43661 0.07152 0.40823
+ 99.000 0.03576 0.50661 0.04655 0.47859
+ 99.900 0.01905 0.58632 0.02634 0.55960
+ 99.990 0.01042 0.64997 0.01530 0.62495
+ 99.999 0.00577 0.70216 0.00901 0.67897
+
+___________________________________________
+2-Sided Confidence Limits For Success Ratio
+___________________________________________
+
+Number of Observations = 200
+Number of successes = 40
+Sample frequency of occurrence = 0.2000000
+
+
+_______________________________________________________________________
+Confidence Lower CP Upper CP Lower JP Upper JP
+ Value (%) Limit Limit Limit Limit
+_______________________________________________________________________
+ 50.000 0.17949 0.22259 0.18190 0.22001
+ 75.000 0.16701 0.23693 0.16934 0.23429
+ 90.000 0.15455 0.25225 0.15681 0.24956
+ 95.000 0.14689 0.26223 0.14910 0.25951
+ 99.000 0.13257 0.28218 0.13468 0.27940
+ 99.900 0.11703 0.30601 0.11902 0.30318
+ 99.990 0.10489 0.32652 0.10677 0.32366
+ 99.999 0.09492 0.34485 0.09670 0.34197
+
+___________________________________________
+2-Sided Confidence Limits For Success Ratio
+___________________________________________
+
+Number of Observations = 2000
+Number of successes = 400
+Sample frequency of occurrence = 0.2000000
+
+
+_______________________________________________________________________
+Confidence Lower CP Upper CP Lower JP Upper JP
+ Value (%) Limit Limit Limit Limit
+_______________________________________________________________________
+ 50.000 0.19382 0.20638 0.19406 0.20613
+ 75.000 0.18965 0.21072 0.18990 0.21047
+ 90.000 0.18537 0.21528 0.18561 0.21503
+ 95.000 0.18267 0.21821 0.18291 0.21796
+ 99.000 0.17745 0.22400 0.17769 0.22374
+ 99.900 0.17150 0.23079 0.17173 0.23053
+ 99.990 0.16658 0.23657 0.16681 0.23631
+ 99.999 0.16233 0.24169 0.16256 0.24143
+
+*/
+
+
+
=======================================
--- /dev/null
+++ /trunk/libs/math/example/binomial_example3.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,91 @@
+// Copyright Paul A. 2006
+// Copyright John Maddock 2006
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Simple example of computing probabilities for a binomial random
variable.
+// Replication of source nag_binomial_dist (g01bjc).
+
+// Shows how to replace NAG C library calls by Boost Math Toolkit C++
calls.
+
+#ifdef _MSC_VER
+# pragma warning(disable: 4512) // assignment operator could not be
generated.
+# pragma warning(disable: 4510) // default constructor could not be
generated.
+# pragma warning(disable: 4610) // can never be instantiated - user
defined constructor required.
+#endif
+
+#include <iostream>
+using std::cout;
+using std::endl;
+using std::ios;
+using std::showpoint;
+#include <iomanip>
+using std::fixed;
+using std::setw;
+#include <boost/math/distributions/binomial.hpp>
+
+int main()
+{
+ cout << "Example 3 of using the binomial distribution to replicate a NAG
library call." << endl;
+ using boost::math::binomial_distribution;
+
+ // This replicates the computation of the examples of using
nag-binomial_dist
+ // using g01bjc in section g01 Sample Calculations on Statistical Data.
+ // http://www.nag.co.uk/numeric/cl/manual/pdf/G01/g01bjc.pdf
+ // Program results section 8.3 page 3.g01bjc.3
+ //8.2. Program Data
+ //g01bjc Example Program Data
+ //4 0.50 2 : n, p, k
+ //19 0.44 13
+ //100 0.75 67
+ //2000 0.33 700
+ //8.3. Program Results
+ //g01bjc Example Program Results
+ //n p k plek pgtk peqk
+ //4 0.500 2 0.68750 0.31250 0.37500
+ //19 0.440 13 0.99138 0.00862 0.01939
+ //100 0.750 67 0.04460 0.95540 0.01700
+ //2000 0.330 700 0.97251 0.02749 0.00312
+
+ cout.setf(ios::showpoint); // Trailing zeros to show significant decimal
digits.
+ cout.precision(5); // Should be able to calculate this?
+ cout << fixed;
+ // Binomial distribution.
+
+ // Note that cdf(dist, k) is equivalent to NAG library plek
probability of <= k
+ // cdf(complement(dist, k)) is equivalent to NAG library pgtk
probability of > k
+ // pdf(dist, k) is equivalent to NAG library peqk
probability of == k
+
+ cout << " n p k plek pgtk peqk " << endl;
+ binomial_distribution<>my_dist(4, 0.5);
+ cout << setw(4) << (int)my_dist.trials() << " " <<
my_dist.success_fraction() << " "<< 2 << " "
+ << cdf(my_dist, 2) << " " << cdf(complement(my_dist, 2)) << " " <<
pdf(my_dist, 2) << endl;
+ binomial_distribution<>two(19, 0.440);
+ cout << setw(4) << (int)two.trials() << " " << two.success_fraction()
<< " " << 13 << " "
+ << cdf(two, 13) << " " << cdf(complement(two, 13)) << " " <<
pdf(two, 13) << endl;
+ binomial_distribution<>three(100, 0.750);
+ cout << setw(4) << (int)three.trials() << " " <<
three.success_fraction() << " " << 67 << " "
+ << cdf(three, 67) << " " << cdf(complement(three, 67)) << " " <<
pdf(three, 67) << endl;
+ binomial_distribution<>four(2000, 0.330);
+ cout << setw(4) << (int)four.trials() << " " <<
four.success_fraction() << " " << 700 << " "
+ << cdf(four, 700) << " " << cdf(complement(four, 700)) << " " <<
pdf(four, 700) << endl;
+
+ return 0;
+} // int main()
+
+/*
+
+Autorun
"i:\boost-sandbox\math_toolkit\libs\math\test\msvc80\debug\binomial_example3.exe"
+Example 3 of using the binomial distribution.
+ n p k plek pgtk peqk
+ 4 0.50000 2 0.68750 0.31250 0.37500
+ 19 0.44000 13 0.99138 0.00862 0.01939
+ 100 0.75000 67 0.04460 0.95540 0.01700
+2000 0.33000 700 0.97251 0.02749 0.00312
+
+
+
+ */
+
=======================================
--- /dev/null
+++ /trunk/libs/math/example/binomial_example_nag.cpp Wed Jul 29 07:28:12
2009
@@ -0,0 +1,91 @@
+// Copyright Paul A. 2007
+// Copyright John Maddock 2007
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Simple example of computing probabilities for a binomial random
variable.
+// Replication of source nag_binomial_dist (g01bjc).
+
+// Shows how to replace NAG C library calls by Boost Math Toolkit C++
calls.
+// Note that the default policy does not replicate the way that NAG
+// library calls handle 'bad' arguments, but you can define policies that
do,
+// as well as other policies that may suit your application even better.
+// See the examples of changing default policies for details.
+
+#include <boost/math/distributions/binomial.hpp>
+
+#include <iostream>
+ using std::cout; using std::endl; using std::ios; using std::showpoint;
+#include <iomanip>
+ using std::fixed; using std::setw;
+
+int main()
+{
+ cout << "Using the binomial distribution to replicate a NAG library
call." << endl;
+ using boost::math::binomial_distribution;
+
+ // This replicates the computation of the examples of using
nag-binomial_dist
+ // using g01bjc in section g01 Somple Calculations on Statistical Data.
+ // http://www.nag.co.uk/numeric/cl/manual/pdf/G01/g01bjc.pdf
+ // Program results section 8.3 page 3.g01bjc.3
+ //8.2. Program Data
+ //g01bjc Example Program Data
+ //4 0.50 2 : n, p, k
+ //19 0.44 13
+ //100 0.75 67
+ //2000 0.33 700
+ //8.3. Program Results
+ //g01bjc Example Program Results
+ //n p k plek pgtk peqk
+ //4 0.500 2 0.68750 0.31250 0.37500
+ //19 0.440 13 0.99138 0.00862 0.01939
+ //100 0.750 67 0.04460 0.95540 0.01700
+ //2000 0.330 700 0.97251 0.02749 0.00312
+
+ cout.setf(ios::showpoint); // Trailing zeros to show significant decimal
digits.
+ cout.precision(5); // Might calculate this from trials in distribution?
+ cout << fixed;
+ // Binomial distribution.
+
+ // Note that cdf(dist, k) is equivalent to NAG library plek
probability of <= k
+ // cdf(complement(dist, k)) is equivalent to NAG library pgtk
probability of > k
+ // pdf(dist, k) is equivalent to NAG library peqk
probability of == k
+
+ cout << " n p k plek pgtk peqk " << endl;
+ binomial_distribution<>my_dist(4, 0.5);
+ cout << setw(4) << (int)my_dist.trials() << " " <<
my_dist.success_fraction()
+ << " " << 2 << " " << cdf(my_dist, 2) << " "
+ << cdf(complement(my_dist, 2)) << " " << pdf(my_dist, 2) << endl;
+
+ binomial_distribution<>two(19, 0.440);
+ cout << setw(4) << (int)two.trials() << " " << two.success_fraction()
+ << " " << 13 << " " << cdf(two, 13) << " "
+ << cdf(complement(two, 13)) << " " << pdf(two, 13) << endl;
+
+ binomial_distribution<>three(100, 0.750);
+ cout << setw(4) << (int)three.trials() << " " <<
three.success_fraction()
+ << " " << 67 << " " << cdf(three, 67) << " " <<
cdf(complement(three, 67))
+ << " " << pdf(three, 67) << endl;
+ binomial_distribution<>four(2000, 0.330);
+ cout << setw(4) << (int)four.trials() << " " <<
four.success_fraction()
+ << " " << 700 << " "
+ << cdf(four, 700) << " " << cdf(complement(four, 700))
+ << " " << pdf(four, 700) << endl;
+
+ return 0;
+} // int main()
+
+/*
+
+Example of using the binomial distribution to replicate a NAG library call.
+ n p k plek pgtk peqk
+ 4 0.50000 2 0.68750 0.31250 0.37500
+ 19 0.44000 13 0.99138 0.00862 0.01939
+ 100 0.75000 67 0.04460 0.95540 0.01700
+2000 0.33000 700 0.97251 0.02749 0.00312
+
+
+ */
+
=======================================
--- /dev/null
+++ /trunk/libs/math/example/binomial_quiz_example.cpp Wed Jul 29 07:28:12
2009
@@ -0,0 +1,520 @@
+// Copyright Paul A. Bristow 2007
+// Copyright John Maddock 2006
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// binomial_examples_quiz.cpp
+
+// Simple example of computing probabilities and quantiles for a binomial
random variable
+// representing the correct guesses on a multiple-choice test.
+
+// source http://www.stat.wvu.edu/SRS/Modules/Binomial/test.html
+
+//[binomial_quiz_example1
+/*`
+A multiple choice test has four possible answers to each of 16 questions.
+A student guesses the answer to each question,
+so the probability of getting a correct answer on any given question is
+one in four, a quarter, 1/4, 25% or fraction 0.25.
+The conditions of the binomial experiment are assumed to be met:
+n = 16 questions constitute the trials;
+each question results in one of two possible outcomes (correct or
incorrect);
+the probability of being correct is 0.25 and is constant if no knowledge
about the subject is assumed;
+the questions are answered independently if the student's answer to a
question
+in no way influences his/her answer to another question.
+
+First, we need to be able to use the binomial distribution constructor
+(and some std input/output, of course).
+*/
+
+#include <boost/math/distributions/binomial.hpp>
+ using boost::math::binomial;
+
+#include <iostream>
+ using std::cout; using std::endl;
+ using std::ios; using std::flush; using std::left; using std::right;
using std::fixed;
+#include <iomanip>
+ using std::setw; using std::setprecision;
+//][/binomial_quiz_example1]
+
+int main()
+{
+ try
+ {
+ cout << "Binomial distribution example - guessing in a quiz." << endl;
+//[binomial_quiz_example2
+/*`
+The number of correct answers, X, is distributed as a binomial random
variable
+with binomial distribution parameters: questions n = 16 and success
fraction probability p = 0.25.
+So we construct a binomial distribution:
+*/
+ int questions = 16; // All the questions in the quiz.
+ int answers = 4; // Possible answers to each question.
+ double success_fraction = (double)answers / (double)questions; // If a
random guess.
+ // Caution: = answers / questions would be zero (because they are
integers)!
+ binomial quiz(questions, success_fraction);
+/*`
+and display the distribution parameters we used thus:
+*/
+ cout << "In a quiz with " << quiz.trials()
+ << " questions and with a probability of guessing right of "
+ << quiz.success_fraction() * 100 << " %"
+ << " or 1 in " << static_cast<int>(1. / quiz.success_fraction()) <<
endl;
+/*`
+Show a few probabilities of just guessing:
+*/
+ cout << "Probability of getting none right is " << pdf(quiz, 0) << endl;
// 0.010023
+ cout << "Probability of getting exactly one right is " << pdf(quiz, 1)
<< endl;
+ cout << "Probability of getting exactly two right is " << pdf(quiz, 2)
<< endl;
+ int pass_score = 11;
+ cout << "Probability of getting exactly " << pass_score << " answers
right by chance is "
+ << pdf(quiz, questions) << endl;
+/*`
+[pre
+Probability of getting none right is 0.0100226
+Probability of getting exactly one right is 0.0534538
+Probability of getting exactly two right is 0.133635
+Probability of getting exactly 11 answers right by chance is 2.32831e-010
+]
+These don't give any encouragement to guessers!
+
+We can tabulate the 'getting exactly right' ( == ) probabilities thus:
+*/
+ cout << "\n" "Guessed Probability" << right << endl;
+ for (int successes = 0; successes <= questions; successes++)
+ {
+ double probability = pdf(quiz, successes);
+ cout << setw(2) << successes << " " << probability << endl;
+ }
+ cout << endl;
+/*`
+[pre
+Guessed Probability
+ 0 0.0100226
+ 1 0.0534538
+ 2 0.133635
+ 3 0.207876
+ 4 0.225199
+ 5 0.180159
+ 6 0.110097
+ 7 0.0524273
+ 8 0.0196602
+ 9 0.00582526
+10 0.00135923
+11 0.000247132
+12 3.43239e-005
+13 3.5204e-006
+14 2.51457e-007
+15 1.11759e-008
+16 2.32831e-010
+]
+Then we can add the probabilities of some 'exactly right' like this:
+*/
+ cout << "Probability of getting none or one right is " << pdf(quiz, 0) +
pdf(quiz, 1) << endl;
+
+/*`
+[pre
+Probability of getting none or one right is 0.0634764
+]
+But if more than a couple of scores are involved, it is more convenient
(and may be more accurate)
+to use the Cumulative Distribution Function (cdf) instead:
+*/
+ cout << "Probability of getting none or one right is " << cdf(quiz, 1)
<< endl;
+/*`
+[pre
+Probability of getting none or one right is 0.0634764
+]
+Since the cdf is inclusive, we can get the probability of getting up to 10
right ( <= )
+*/
+ cout << "Probability of getting <= 10 right (to fail) is " << cdf(quiz,
10) << endl;
+/*`
+[pre
+Probability of getting <= 10 right (to fail) is 0.999715
+]
+To get the probability of getting 11 or more right (to pass),
+it is tempting to use ``1 - cdf(quiz, 10)`` to get the probability of > 10
+*/
+ cout << "Probability of getting > 10 right (to pass) is " << 1 -
cdf(quiz, 10) << endl;
+/*`
+[pre
+Probability of getting > 10 right (to pass) is 0.000285239
+]
+But this should be resisted in favor of using the complement function.
[link why_complements Why complements?]
+*/
+ cout << "Probability of getting > 10 right (to pass) is " <<
cdf(complement(quiz, 10)) << endl;
+/*`
+[pre
+Probability of getting > 10 right (to pass) is 0.000285239
+]
+And we can check that these two, <= 10 and > 10, add up to unity.
+*/
+BOOST_ASSERT((cdf(quiz, 10) + cdf(complement(quiz, 10))) == 1.);
+/*`
+If we want a < rather than a <= test, because the CDF is inclusive, we
must subtract one from the score.
+*/
+ cout << "Probability of getting less than " << pass_score
+ << " (< " << pass_score << ") answers right by guessing is "
+ << cdf(quiz, pass_score -1) << endl;
+/*`
+[pre
+Probability of getting less than 11 (< 11) answers right by guessing is
0.999715
+]
+and similarly to get a >= rather than a > test
+we also need to subtract one from the score (and can again check the sum
is unity).
+This is because if the cdf is /inclusive/,
+then its complement must be /exclusive/ otherwise there would be one
possible
+outcome counted twice!
+*/
+ cout << "Probability of getting at least " << pass_score
+ << "(>= " << pass_score << ") answers right by guessing is "
+ << cdf(complement(quiz, pass_score-1))
+ << ", only 1 in " << 1/cdf(complement(quiz, pass_score-1)) << endl;
+
+ BOOST_ASSERT((cdf(quiz, pass_score -1) + cdf(complement(quiz,
pass_score-1))) == 1);
+
+/*`
+[pre
+Probability of getting at least 11 (>= 11) answers right by guessing is
0.000285239, only 1 in 3505.83
+]
+Finally we can tabulate some probabilities:
+*/
+ cout << "\n" "At most (<=)""\n""Guessed OK Probability" << right <<
endl;
+ for (int score = 0; score <= questions; score++)
+ {
+ cout << setw(2) << score << " " << setprecision(10)
+ << cdf(quiz, score) << endl;
+ }
+ cout << endl;
+/*`
+[pre
+At most (<=)
+Guessed OK Probability
+ 0 0.01002259576
+ 1 0.0634764398
+ 2 0.1971110499
+ 3 0.4049871101
+ 4 0.6301861752
+ 5 0.8103454274
+ 6 0.9204427481
+ 7 0.9728700437
+ 8 0.9925302796
+ 9 0.9983555346
+10 0.9997147608
+11 0.9999618928
+12 0.9999962167
+13 0.9999997371
+14 0.9999999886
+15 0.9999999998
+16 1
+]
+*/
+ cout << "\n" "At least (>)""\n""Guessed OK Probability" << right <<
endl;
+ for (int score = 0; score <= questions; score++)
+ {
+ cout << setw(2) << score << " " << setprecision(10)
+ << cdf(complement(quiz, score)) << endl;
+ }
+/*`
+[pre
+At least (>)
+Guessed OK Probability
+ 0 0.9899774042
+ 1 0.9365235602
+ 2 0.8028889501
+ 3 0.5950128899
+ 4 0.3698138248
+ 5 0.1896545726
+ 6 0.07955725188
+ 7 0.02712995629
+ 8 0.00746972044
+ 9 0.001644465374
+10 0.0002852391917
+11 3.810715862e-005
+12 3.783265129e-006
+13 2.628657967e-007
+14 1.140870154e-008
+15 2.328306437e-010
+16 0
+]
+We now consider the probabilities of *ranges* of correct guesses.
+
+First, calculate the probability of getting a range of guesses right,
+by adding the exact probabilities of each from low ... high.
+*/
+ int low = 3; // Getting at least 3 right.
+ int high = 5; // Getting as most 5 right.
+ double sum = 0.;
+ for (int i = low; i <= high; i++)
+ {
+ sum += pdf(quiz, i);
+ }
+ cout.precision(4);
+ cout << "Probability of getting between "
+ << low << " and " << high << " answers right by guessing is "
+ << sum << endl; // 0.61323
+/*`
+[pre
+Probability of getting between 3 and 5 answers right by guessing is 0.6132
+]
+Or, usually better, we can use the difference of cdfs instead:
+*/
+ cout << "Probability of getting between " << low << " and " << high << "
answers right by guessing is "
+ << cdf(quiz, high) - cdf(quiz, low - 1) << endl; // 0.61323
+/*`
+[pre
+Probability of getting between 3 and 5 answers right by guessing is 0.6132
+]
+And we can also try a few more combinations of high and low choices:
+*/
+ low = 1; high = 6;
+ cout << "Probability of getting between " << low << " and " << high << "
answers right by guessing is "
+ << cdf(quiz, high) - cdf(quiz, low - 1) << endl; // 1 and 6 P= 0.91042
+ low = 1; high = 8;
+ cout << "Probability of getting between " << low << " and " << high << "
answers right by guessing is "
+ << cdf(quiz, high) - cdf(quiz, low - 1) << endl; // 1 <= x 8 P =
0.9825
+ low = 4; high = 4;
+ cout << "Probability of getting between " << low << " and " << high << "
answers right by guessing is "
+ << cdf(quiz, high) - cdf(quiz, low - 1) << endl; // 4 <= x 4 P =
0.22520
+
+/*`
+[pre
+Probability of getting between 1 and 6 answers right by guessing is 0.9104
+Probability of getting between 1 and 8 answers right by guessing is 0.9825
+Probability of getting between 4 and 4 answers right by guessing is 0.2252
+]
+[h4 Using Binomial distribution moments]
+Using moments of the distribution, we can say more about the spread of
results from guessing.
+*/
+ cout << "By guessing, on average, one can expect to get " << mean(quiz)
<< " correct answers." << endl;
+ cout << "Standard deviation is " << standard_deviation(quiz) << endl;
+ cout << "So about 2/3 will lie within 1 standard deviation and get
between "
+ << ceil(mean(quiz) - standard_deviation(quiz)) << " and "
+ << floor(mean(quiz) + standard_deviation(quiz)) << " correct." << endl;
+ cout << "Mode (the most frequent) is " << mode(quiz) << endl;
+ cout << "Skewness is " << skewness(quiz) << endl;
+
+/*`
+[pre
+By guessing, on average, one can expect to get 4 correct answers.
+Standard deviation is 1.732
+So about 2/3 will lie within 1 standard deviation and get between 3 and 5
correct.
+Mode (the most frequent) is 4
+Skewness is 0.2887
+]
+[h4 Quantiles]
+The quantiles (percentiles or percentage points) for a few probability
levels:
+*/
+ cout << "Quartiles " << quantile(quiz, 0.25) << " to "
+ << quantile(complement(quiz, 0.25)) << endl; // Quartiles
+ cout << "1 standard deviation " << quantile(quiz, 0.33) << " to "
+ << quantile(quiz, 0.67) << endl; // 1 sd
+ cout << "Deciles " << quantile(quiz, 0.1) << " to "
+ << quantile(complement(quiz, 0.1))<< endl; // Deciles
+ cout << "5 to 95% " << quantile(quiz, 0.05) << " to "
+ << quantile(complement(quiz, 0.05))<< endl; // 5 to 95%
+ cout << "2.5 to 97.5% " << quantile(quiz, 0.025) << " to "
+ << quantile(complement(quiz, 0.025)) << endl; // 2.5 to 97.5%
+ cout << "2 to 98% " << quantile(quiz, 0.02) << " to "
+ << quantile(complement(quiz, 0.02)) << endl; // 2 to 98%
+
+ cout << "If guessing then percentiles 1 to 99% will get " <<
quantile(quiz, 0.01)
+ << " to " << quantile(complement(quiz, 0.01)) << " right." << endl;
+/*`
+Notice that these output integral values because the default policy is
`integer_round_outwards`.
+[pre
+Quartiles 2 to 5
+1 standard deviation 2 to 5
+Deciles 1 to 6
+5 to 95% 0 to 7
+2.5 to 97.5% 0 to 8
+2 to 98% 0 to 8
+]
+*/
+
+//] [/binomial_quiz_example2]
+
+//[discrete_quantile_real
+/*`
+Quantiles values are controlled by the
+[link math_toolkit.policy.pol_ref.discrete_quant_ref discrete quantile
policy]
+chosen.
+The default is `integer_round_outwards`,
+so the lower quantile is rounded down, and the upper quantile is rounded
up.
+
+But we might believe that the real values tell us a little more - see
+[link math_toolkit.policy.pol_tutorial.understand_dis_quant Understanding
Discrete Quantile Policy].
+
+We could control the policy for *all* distributions by
+
+ #define BOOST_MATH_DISCRETE_QUANTILE_POLICY real
+
+ at the head of the program would make this policy apply
+to this *one, and only*, translation unit.
+
+Or we can now create a (typedef for) policy that has discrete quantiles
real
+(here avoiding any 'using namespaces ...' statements):
+*/
+ using boost::math::policies::policy;
+ using boost::math::policies::discrete_quantile;
+ using boost::math::policies::real;
+ using boost::math::policies::integer_round_outwards; // Default.
+ typedef boost::math::policies::policy<discrete_quantile<real> >
real_quantile_policy;
+/*`
+Add a custom binomial distribution called ``real_quantile_binomial`` that
uses ``real_quantile_policy``
+*/
+ using boost::math::binomial_distribution;
+ typedef binomial_distribution<double, real_quantile_policy>
real_quantile_binomial;
+/*`
+Construct an object of this custom distribution:
+*/
+ real_quantile_binomial quiz_real(questions, success_fraction);
+/*`
+And use this to show some quantiles - that now have real rather than
integer values.
+*/
+ cout << "Quartiles " << quantile(quiz, 0.25) << " to "
+ << quantile(complement(quiz_real, 0.25)) << endl; // Quartiles 2 to
4.6212
+ cout << "1 standard deviation " << quantile(quiz_real, 0.33) << " to "
+ << quantile(quiz_real, 0.67) << endl; // 1 sd 2.6654 4.194
+ cout << "Deciles " << quantile(quiz_real, 0.1) << " to "
+ << quantile(complement(quiz_real, 0.1))<< endl; // Deciles 1.3487
5.7583
+ cout << "5 to 95% " << quantile(quiz_real, 0.05) << " to "
+ << quantile(complement(quiz_real, 0.05))<< endl; // 5 to 95% 0.83739
6.4559
+ cout << "2.5 to 97.5% " << quantile(quiz_real, 0.025) << " to "
+ << quantile(complement(quiz_real, 0.025)) << endl; // 2.5 to 97.5%
0.42806 7.0688
+ cout << "2 to 98% " << quantile(quiz_real, 0.02) << " to "
+ << quantile(complement(quiz_real, 0.02)) << endl; // 2 to 98% 0.31311
7.7880
+
+ cout << "If guessing, then percentiles 1 to 99% will get " <<
quantile(quiz_real, 0.01)
+ << " to " << quantile(complement(quiz_real, 0.01)) << " right." <<
endl;
+/*`
+[pre
+Real Quantiles
+Quartiles 2 to 4.621
+1 standard deviation 2.665 to 4.194
+Deciles 1.349 to 5.758
+5 to 95% 0.8374 to 6.456
+2.5 to 97.5% 0.4281 to 7.069
+2 to 98% 0.3131 to 7.252
+If guessing then percentiles 1 to 99% will get 0 to 7.788 right.
+]
+*/
+
+//] [/discrete_quantile_real]
+ }
+ catch(const std::exception& e)
+ { // Always useful to include try & catch blocks because
+ // default policies are to throw exceptions on arguments that cause
+ // errors like underflow, overflow.
+ // Lacking try & catch blocks, the program will abort without a
message below,
+ // which may give some helpful clues as to the cause of the exception.
+ std::cout <<
+ "\n""Message from thrown exception was:\n " << e.what() <<
std::endl;
+ }
+ return 0;
+} // int main()
+
+
+
+/*
+
+Output is:
+
+Binomial distribution example - guessing in a quiz.
+In a quiz with 16 questions and with a probability of guessing right of
25 % or 1 in 4
+Probability of getting none right is 0.0100226
+Probability of getting exactly one right is 0.0534538
+Probability of getting exactly two right is 0.133635
+Probability of getting exactly 11 answers right by chance is 2.32831e-010
+Guessed Probability
+ 0 0.0100226
+ 1 0.0534538
+ 2 0.133635
+ 3 0.207876
+ 4 0.225199
+ 5 0.180159
+ 6 0.110097
+ 7 0.0524273
+ 8 0.0196602
+ 9 0.00582526
+10 0.00135923
+11 0.000247132
+12 3.43239e-005
+13 3.5204e-006
+14 2.51457e-007
+15 1.11759e-008
+16 2.32831e-010
+Probability of getting none or one right is 0.0634764
+Probability of getting none or one right is 0.0634764
+Probability of getting <= 10 right (to fail) is 0.999715
+Probability of getting > 10 right (to pass) is 0.000285239
+Probability of getting > 10 right (to pass) is 0.000285239
+Probability of getting less than 11 (< 11) answers right by guessing is
0.999715
+Probability of getting at least 11(>= 11) answers right by guessing is
0.000285239, only 1 in 3505.83
+At most (<=)
+Guessed OK Probability
+ 0 0.01002259576
+ 1 0.0634764398
+ 2 0.1971110499
+ 3 0.4049871101
+ 4 0.6301861752
+ 5 0.8103454274
+ 6 0.9204427481
+ 7 0.9728700437
+ 8 0.9925302796
+ 9 0.9983555346
+10 0.9997147608
+11 0.9999618928
+12 0.9999962167
+13 0.9999997371
+14 0.9999999886
+15 0.9999999998
+16 1
+At least (>=)
+Guessed OK Probability
+ 0 0.9899774042
+ 1 0.9365235602
+ 2 0.8028889501
+ 3 0.5950128899
+ 4 0.3698138248
+ 5 0.1896545726
+ 6 0.07955725188
+ 7 0.02712995629
+ 8 0.00746972044
+ 9 0.001644465374
+10 0.0002852391917
+11 3.810715862e-005
+12 3.783265129e-006
+13 2.628657967e-007
+14 1.140870154e-008
+15 2.328306437e-010
+16 0
+Probability of getting between 3 and 5 answers right by guessing is 0.6132
+Probability of getting between 3 and 5 answers right by guessing is 0.6132
+Probability of getting between 1 and 6 answers right by guessing is 0.9104
+Probability of getting between 1 and 8 answers right by guessing is 0.9825
+Probability of getting between 4 and 4 answers right by guessing is 0.2252
+By guessing, on average, one can expect to get 4 correct answers.
+Standard deviation is 1.732
+So about 2/3 will lie within 1 standard deviation and get between 3 and 5
correct.
+Mode (the most frequent) is 4
+Skewness is 0.2887
+Quartiles 2 to 5
+1 standard deviation 2 to 5
+Deciles 1 to 6
+5 to 95% 0 to 7
+2.5 to 97.5% 0 to 8
+2 to 98% 0 to 8
+If guessing then percentiles 1 to 99% will get 0 to 8 right.
+Quartiles 2 to 4.621
+1 standard deviation 2.665 to 4.194
+Deciles 1.349 to 5.758
+5 to 95% 0.8374 to 6.456
+2.5 to 97.5% 0.4281 to 7.069
+2 to 98% 0.3131 to 7.252
+If guessing, then percentiles 1 to 99% will get 0 to 7.788 right.
+
+*/
+
=======================================
--- /dev/null
+++ /trunk/libs/math/example/binomial_sample_sizes.cpp Wed Jul 29 07:28:12
2009
@@ -0,0 +1,76 @@
+// (C) Copyright John Maddock 2006
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifdef _MSC_VER
+# pragma warning(disable: 4512) // assignment operator could not be
generated.
+# pragma warning(disable: 4510) // default constructor could not be
generated.
+# pragma warning(disable: 4610) // can never be instantiated - user
defined constructor required.
+#endif
+
+#include <iostream>
+#include <iomanip>
+#include <boost/math/distributions/binomial.hpp>
+
+void find_max_sample_size(double p, unsigned successes)
+{
+ //
+ // p = success ratio.
+ // successes = Total number of observed successes.
+ //
+ // Calculate how many trials we can have to ensure the
+ // maximum number of successes does not exceed "successes".
+ // A typical use would be failure analysis, where you want
+ // zero or fewer "successes" with some probability.
+ //
+ using namespace std;
+ using namespace boost::math;
+
+ // Print out general info:
+ cout <<
+ "________________________\n"
+ "Maximum Number of Trials\n"
+ "________________________\n\n";
+ cout << setprecision(7);
+ cout << setw(40) << left << "Success ratio" << "= " << p << "\n";
+ cout << setw(40) << left << "Maximum Number of \"successes\" permitted"
<< "= " << successes << "\n";
+ //
+ // Define a table of confidence intervals:
+ //
+ double alpha[] = { 0.5, 0.25, 0.1, 0.05, 0.01, 0.001, 0.0001, 0.00001 };
+ //
+ // Print table header:
+ //
+ cout << "\n\n"
+ "____________________________\n"
+ "Confidence Max Number\n"
+ " Value (%) Of Trials \n"
+ "____________________________\n";
+ //
+ // Now print out the data for the table rows.
+ //
+ for(unsigned i = 0; i < sizeof(alpha)/sizeof(alpha[0]); ++i)
+ {
+ // Confidence value:
+ cout << fixed << setprecision(3) << setw(10) << right << 100 *
(1-alpha[i]);
+ // calculate trials:
+ double t =
binomial_distribution<>::find_maximum_number_of_trials(successes, p,
alpha[i]);
+ t = floor(t);
+ // Print Trials:
+ cout << fixed << setprecision(0) << setw(15) << right << t << endl;
+ }
+ cout << endl;
+}
+
+int main()
+{
+ find_max_sample_size(1.0/1000, 0);
+ find_max_sample_size(1.0/10000, 0);
+ find_max_sample_size(1.0/100000, 0);
+ find_max_sample_size(1.0/1000000, 0);
+
+ return 0;
+}
+
=======================================
--- /dev/null
+++ /trunk/libs/math/example/c_error_policy_example.cpp Wed Jul 29 07:28:12
2009
@@ -0,0 +1,74 @@
+// C_error_policy_example.cpp
+
+// Copyright Paul A. Bristow 2007.
+// Copyright John Maddock 2007.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Suppose we want a call to tgamma to behave in a C-compatible way
+// and set global ::errno rather than throw an exception.
+
+#include <boost/math/special_functions/gamma.hpp>
+using boost::math::tgamma;
+
+using boost::math::policies::policy;
+// Possible errors
+using boost::math::policies::overflow_error;
+using boost::math::policies::underflow_error;
+using boost::math::policies::domain_error;
+using boost::math::policies::pole_error;
+using boost::math::policies::denorm_error;
+using boost::math::policies::evaluation_error;
+
+using boost::math::policies::errno_on_error;
+//using boost::math::policies::ignore_error;
+
+//using namespace boost::math::policies;
+//using namespace boost::math;
+
+// Define a policy:
+typedef policy<
+ domain_error<errno_on_error>, // 'bad' arguments.
+ pole_error<errno_on_error>, // argument is pole value.
+ overflow_error<errno_on_error>, // argument value causes overflow.
+ evaluation_error<errno_on_error> // evaluation does not converge
and may be inaccurate, or worse.
+ > C_error_policy;
+
+// std
+#include <iostream>
+ using std::cout;
+ using std::endl;
+
+int main()
+{
+ // We can achieve this at the function call site
+ // with the previously defined policy C_error_policy.
+ double t = tgamma(4., C_error_policy());
+ cout << "tgamma(4., C_error_policy() = " << t << endl; // 6
+
+ // Alternatively we could use the function make_policy,
+ // provided for convenience,
+ // and define everything at the call site:
+ t = tgamma(4., make_policy(
+ domain_error<errno_on_error>(),
+ pole_error<errno_on_error>(),
+ overflow_error<errno_on_error>(),
+ evaluation_error<errno_on_error>()
+ ));
+ cout << "tgamma(4., make_policy( ...) = " << t << endl; // 6
+
+ return 0;
+} // int main()
+
+/*
+
+Output
+
+Autorun
"i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\c_error_policy_example.exe"
+tgamma(4., C_error_policy() = 6
+tgamma(4., make_policy( ...) = 6
+
+*/
=======================================
--- /dev/null
+++ /trunk/libs/math/example/chi_square_std_dev_test.cpp Wed Jul 29
07:28:12 2009
@@ -0,0 +1,552 @@
+// Copyright John Maddock 2006, 2007
+// Copyright Paul A. Bristow 2007
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include <iostream>
+using std::cout; using std::endl;
+using std::left; using std::fixed; using std::right; using std::scientific;
+#include <iomanip>
+using std::setw;
+using std::setprecision;
+#include <boost/math/distributions/chi_squared.hpp>
+
+
+void confidence_limits_on_std_deviation(
+ double Sd, // Sample Standard Deviation
+ unsigned N) // Sample size
+{
+ // Calculate confidence intervals for the standard deviation.
+ // For example if we set the confidence limit to
+ // 0.95, we know that if we repeat the sampling
+ // 100 times, then we expect that the true standard deviation
+ // will be between out limits on 95 occations.
+ // Note: this is not the same as saying a 95%
+ // confidence interval means that there is a 95%
+ // probability that the interval contains the true standard deviation.
+ // The interval computed from a given sample either
+ // contains the true standard deviation or it does not.
+ // See http://www.itl.nist.gov/div898/handbook/eda/section3/eda358.htm
+
+ // using namespace boost::math;
+ using boost::math::chi_squared;
+ using boost::math::quantile;
+ using boost::math::complement;
+
+ // Print out general info:
+ cout <<
+ "________________________________________________\n"
+ "2-Sided Confidence Limits For Standard Deviation\n"
+ "________________________________________________\n\n";
+ cout << setprecision(7);
+ cout << setw(40) << left << "Number of Observations" << "= " << N
<< "\n";
+ cout << setw(40) << left << "Standard Deviation" << "= " << Sd << "\n";
+ //
+ // Define a table of significance/risk levels:
+ double alpha[] = { 0.5, 0.25, 0.1, 0.05, 0.01, 0.001, 0.0001, 0.00001 };
+ //
+ // Start by declaring the distribution we'll need:
+ chi_squared dist(N - 1);
+ //
+ // Print table header:
+ //
+ cout << "\n\n"
+ "_____________________________________________\n"
+ "Confidence Lower Upper\n"
+ " Value (%) Limit Limit\n"
+ "_____________________________________________\n";
+ //
+ // Now print out the data for the table rows.
+ for(unsigned i = 0; i < sizeof(alpha)/sizeof(alpha[0]); ++i)
+ {
+ // Confidence value:
+ cout << fixed << setprecision(3) << setw(10) << right << 100 *
(1-alpha[i]);
+ // Calculate limits:
+ double lower_limit = sqrt((N - 1) * Sd * Sd /
quantile(complement(dist, alpha[i] / 2)));
+ double upper_limit = sqrt((N - 1) * Sd * Sd / quantile(dist,
alpha[i] / 2));
+ // Print Limits:
+ cout << fixed << setprecision(5) << setw(15) << right << lower_limit;
+ cout << fixed << setprecision(5) << setw(15) << right << upper_limit
<< endl;
+ }
+ cout << endl;
+} // void confidence_limits_on_std_deviation
+
+void confidence_limits_on_std_deviation_alpha(
+ double Sd, // Sample Standard Deviation
+ double alpha // confidence
+ )
+{ // Calculate confidence intervals for the standard deviation.
+ // for the alpha parameter, for a range number of observations,
+ // from a mere 2 up to a million.
+ // O. L. Davies, Statistical Methods in Research and Production, ISBN 0
05 002437 X,
+ // 4.33 Page 68, Table H, pp 452 459.
+
+ // using namespace std;
+ // using namespace boost::math;
+ using boost::math::chi_squared;
+ using boost::math::quantile;
+ using boost::math::complement;
+
+ // Define a table of numbers of observations:
+ unsigned int obs[] = {2, 3, 4, 5, 6, 7, 8, 9, 10, 15, 20, 30, 40 , 50,
60, 100, 120, 1000, 10000, 50000, 100000, 1000000};
+
+ cout << // Print out heading:
+ "________________________________________________\n"
+ "2-Sided Confidence Limits For Standard Deviation\n"
+ "________________________________________________\n\n";
+ cout << setprecision(7);
+ cout << setw(40) << left << "Confidence level (two-sided) " << "= " <<
alpha << "\n";
+ cout << setw(40) << left << "Standard Deviation" << "= " << Sd << "\n";
+
+ cout << "\n\n" // Print table header:
+ "_____________________________________________\n"
+ "Observations Lower Upper\n"
+ " Limit Limit\n"
+ "_____________________________________________\n";
+ for(unsigned i = 0; i < sizeof(obs)/sizeof(obs[0]); ++i)
+ {
+ unsigned int N = obs[i]; // Observations
+ // Start by declaring the distribution with the appropriate :
+ chi_squared dist(N - 1);
+
+ // Now print out the data for the table row.
+ cout << fixed << setprecision(3) << setw(10) << right << N;
+ // Calculate limits: (alpha /2 because it is a two-sided (upper and
lower limit) test.
+ double lower_limit = sqrt((N - 1) * Sd * Sd /
quantile(complement(dist, alpha / 2)));
+ double upper_limit = sqrt((N - 1) * Sd * Sd / quantile(dist, alpha /
2));
+ // Print Limits:
+ cout << fixed << setprecision(4) << setw(15) << right << lower_limit;
+ cout << fixed << setprecision(4) << setw(15) << right << upper_limit
<< endl;
+ }
+ cout << endl;
+}// void confidence_limits_on_std_deviation_alpha
+
+void chi_squared_test(
+ double Sd, // Sample std deviation
+ double D, // True std deviation
+ unsigned N, // Sample size
+ double alpha) // Significance level
+{
+ //
+ // A Chi Squared test applied to a single set of data.
+ // We are testing the null hypothesis that the true
+ // standard deviation of the sample is D, and that any variation is down
+ // to chance. We can also test the alternative hypothesis
+ // that any difference is not down to chance.
+ // See http://www.itl.nist.gov/div898/handbook/eda/section3/eda358.htm
+ //
+ // using namespace boost::math;
+ using boost::math::chi_squared;
+ using boost::math::quantile;
+ using boost::math::complement;
+ using boost::math::cdf;
+
+ // Print header:
+ cout <<
+ "______________________________________________\n"
+ "Chi Squared test for sample standard deviation\n"
+ "______________________________________________\n\n";
+ cout << setprecision(5);
+ cout << setw(55) << left << "Number of Observations" << "= " << N
<< "\n";
+ cout << setw(55) << left << "Sample Standard Deviation" << "= " << Sd
<< "\n";
+ cout << setw(55) << left << "Expected True Standard Deviation" << "= "
<< D << "\n\n";
+ //
+ // Now we can calculate and output some stats:
+ //
+ // test-statistic:
+ double t_stat = (N - 1) * (Sd / D) * (Sd / D);
+ cout << setw(55) << left << "Test Statistic" << "= " << t_stat << "\n";
+ //
+ // Finally define our distribution, and get the probability:
+ //
+ chi_squared dist(N - 1);
+ double p = cdf(dist, t_stat);
+ cout << setw(55) << left << "CDF of test statistic: " << "= "
+ << setprecision(3) << scientific << p << "\n";
+ double ucv = quantile(complement(dist, alpha));
+ double ucv2 = quantile(complement(dist, alpha / 2));
+ double lcv = quantile(dist, alpha);
+ double lcv2 = quantile(dist, alpha / 2);
+ cout << setw(55) << left << "Upper Critical Value at alpha: " << "= "
+ << setprecision(3) << scientific << ucv << "\n";
+ cout << setw(55) << left << "Upper Critical Value at alpha/2: " << "= "
+ << setprecision(3) << scientific << ucv2 << "\n";
+ cout << setw(55) << left << "Lower Critical Value at alpha: " << "= "
+ << setprecision(3) << scientific << lcv << "\n";
+ cout << setw(55) << left << "Lower Critical Value at alpha/2: " << "= "
+ << setprecision(3) << scientific << lcv2 << "\n\n";
+ //
+ // Finally print out results of alternative hypothesis:
+ //
+ cout << setw(55) << left <<
+ "Results for Alternative Hypothesis and alpha" << "= "
+ << setprecision(4) << fixed << alpha << "\n\n";
+ cout << "Alternative Hypothesis Conclusion\n";
+ cout << "Standard Deviation != " << setprecision(3) << fixed << D
<< " ";
+ if((ucv2 < t_stat) || (lcv2 > t_stat))
+ cout << "NOT REJECTED\n";
+ else
+ cout << "REJECTED\n";
+ cout << "Standard Deviation < " << setprecision(3) << fixed << D
<< " ";
+ if(lcv > t_stat)
+ cout << "NOT REJECTED\n";
+ else
+ cout << "REJECTED\n";
+ cout << "Standard Deviation > " << setprecision(3) << fixed << D
<< " ";
+ if(ucv < t_stat)
+ cout << "NOT REJECTED\n";
+ else
+ cout << "REJECTED\n";
+ cout << endl << endl;
+} // void chi_squared_test
+
+void chi_squared_sample_sized(
+ double diff, // difference from variance to detect
+ double variance) // true variance
+{
+ using namespace std;
+ // using boost::math;
+ using boost::math::chi_squared;
+ using boost::math::quantile;
+ using boost::math::complement;
+ using boost::math::cdf;
+
+ try
+ {
+ cout << // Print out general info:
+ "_____________________________________________________________\n"
+ "Estimated sample sizes required for various confidence levels\n"
+ "_____________________________________________________________\n\n";
+ cout << setprecision(5);
+ cout << setw(40) << left << "True Variance" << "= " << variance
<< "\n";
+ cout << setw(40) << left << "Difference to detect" << "= " << diff
<< "\n";
+ //
+ // Define a table of significance levels:
+ //
+ double alpha[] = { 0.5, 0.33333333333333333333333, 0.25, 0.1, 0.05,
0.01, 0.001, 0.0001, 0.00001 };
+ //
+ // Print table header:
+ //
+ cout << "\n\n"
+ "_______________________________________________________________\n"
+ "Confidence Estimated Estimated\n"
+ " Value (%) Sample Size Sample Size\n"
+ " (lower one- (upper one-\n"
+ " sided test) sided test)\n"
+ "_______________________________________________________________\n";
+ //
+ // Now print out the data for the table rows.
+ //
+ for(unsigned i = 0; i < sizeof(alpha)/sizeof(alpha[0]); ++i)
+ {
+ // Confidence value:
+ cout << fixed << setprecision(3) << setw(10) << right << 100 *
(1-alpha[i]);
+ // calculate df for a lower single-sided test:
+ double df = chi_squared::find_degrees_of_freedom(
+ -diff, alpha[i], alpha[i], variance);
+ // convert to sample size:
+ double size = ceil(df) + 1;
+ // Print size:
+ cout << fixed << setprecision(0) << setw(16) << right << size;
+ // calculate df for an upper single-sided test:
+ df = chi_squared::find_degrees_of_freedom(
+ diff, alpha[i], alpha[i], variance);
+ // convert to sample size:
+ size = ceil(df) + 1;
+ // Print size:
+ cout << fixed << setprecision(0) << setw(16) << right << size <<
endl;
+ }
+ cout << endl;
+ }
+ catch(const std::exception& e)
+ { // Always useful to include try & catch blocks because default policies
+ // are to throw exceptions on arguments that cause errors like
underflow, overflow.
+ // Lacking try & catch blocks, the program will abort without a
message below,
+ // which may give some helpful clues as to the cause of the exception.
+ std::cout <<
+ "\n""Message from thrown exception was:\n " << e.what() <<
std::endl;
+ }
+} // chi_squared_sample_sized
+
+int main()
+{
+ // Run tests for Gear data
+ // see http://www.itl.nist.gov/div898/handbook/eda/section3/eda3581.htm
+ // Tests measurements of gear diameter.
+ //
+ confidence_limits_on_std_deviation(0.6278908E-02, 100);
+ chi_squared_test(0.6278908E-02, 0.1, 100, 0.05);
+ chi_squared_sample_sized(0.1 - 0.6278908E-02, 0.1);
+ //
+ // Run tests for silicon wafer fabrication data.
+ // see http://www.itl.nist.gov/div898/handbook/prc/section2/prc23.htm
+ // A supplier of 100 ohm.cm silicon wafers claims that his fabrication
+ // process can produce wafers with sufficient consistency so that the
+ // standard deviation of resistivity for the lot does not exceed
+ // 10 ohm.cm. A sample of N = 10 wafers taken from the lot has a
+ // standard deviation of 13.97 ohm.cm
+ //
+ confidence_limits_on_std_deviation(13.97, 10);
+ chi_squared_test(13.97, 10.0, 10, 0.05);
+ chi_squared_sample_sized(13.97 * 13.97 - 100, 100);
+ chi_squared_sample_sized(55, 100);
+ chi_squared_sample_sized(1, 100);
+
+ // List confidence interval multipliers for standard deviation
+ // for a range of numbers of observations from 2 to a million,
+ // and for a few alpha values, 0.1, 0.05, 0.01 for condfidences 90, 95,
99 %
+ confidence_limits_on_std_deviation_alpha(1., 0.1);
+ confidence_limits_on_std_deviation_alpha(1., 0.05);
+ confidence_limits_on_std_deviation_alpha(1., 0.01);
+
+ return 0;
+}
+
+/*
+
+________________________________________________
+2-Sided Confidence Limits For Standard Deviation
+________________________________________________
+Number of Observations = 100
+Standard Deviation = 0.006278908
+_____________________________________________
+Confidence Lower Upper
+ Value (%) Limit Limit
+_____________________________________________
+ 50.000 0.00601 0.00662
+ 75.000 0.00582 0.00685
+ 90.000 0.00563 0.00712
+ 95.000 0.00551 0.00729
+ 99.000 0.00530 0.00766
+ 99.900 0.00507 0.00812
+ 99.990 0.00489 0.00855
+ 99.999 0.00474 0.00895
+______________________________________________
+Chi Squared test for sample standard deviation
+______________________________________________
+Number of Observations = 100
+Sample Standard Deviation = 0.00628
+Expected True Standard Deviation = 0.10000
+Test Statistic = 0.39030
+CDF of test statistic: = 1.438e-099
+Upper Critical Value at alpha: = 1.232e+002
+Upper Critical Value at alpha/2: = 1.284e+002
+Lower Critical Value at alpha: = 7.705e+001
+Lower Critical Value at alpha/2: = 7.336e+001
+Results for Alternative Hypothesis and alpha = 0.0500
+Alternative Hypothesis Conclusion
+Standard Deviation != 0.100 NOT REJECTED
+Standard Deviation < 0.100 NOT REJECTED
+Standard Deviation > 0.100 REJECTED
+_____________________________________________________________
+Estimated sample sizes required for various confidence levels
+_____________________________________________________________
+True Variance = 0.10000
+Difference to detect = 0.09372
+_______________________________________________________________
+Confidence Estimated Estimated
+ Value (%) Sample Size Sample Size
+ (lower one- (upper one-
+ sided test) sided test)
+_______________________________________________________________
+ 50.000 2 2
+ 66.667 2 5
+ 75.000 2 10
+ 90.000 4 32
+ 95.000 5 52
+ 99.000 8 102
+ 99.900 13 178
+ 99.990 18 257
+ 99.999 23 337
+________________________________________________
+2-Sided Confidence Limits For Standard Deviation
+________________________________________________
+Number of Observations = 10
+Standard Deviation = 13.9700000
+_____________________________________________
+Confidence Lower Upper
+ Value (%) Limit Limit
+_____________________________________________
+ 50.000 12.41880 17.25579
+ 75.000 11.23084 19.74131
+ 90.000 10.18898 22.98341
+ 95.000 9.60906 25.50377
+ 99.000 8.62898 31.81825
+ 99.900 7.69466 42.51593
+ 99.990 7.04085 55.93352
+ 99.999 6.54517 73.00132
+______________________________________________
+Chi Squared test for sample standard deviation
+______________________________________________
+Number of Observations = 10
+Sample Standard Deviation = 13.97000
+Expected True Standard Deviation = 10.00000
+Test Statistic = 17.56448
+CDF of test statistic: = 9.594e-001
+Upper Critical Value at alpha: = 1.692e+001
+Upper Critical Value at alpha/2: = 1.902e+001
+Lower Critical Value at alpha: = 3.325e+000
+Lower Critical Value at alpha/2: = 2.700e+000
+Results for Alternative Hypothesis and alpha = 0.0500
+Alternative Hypothesis Conclusion
+Standard Deviation != 10.000 REJECTED
+Standard Deviation < 10.000 REJECTED
+Standard Deviation > 10.000 NOT REJECTED
+_____________________________________________________________
+Estimated sample sizes required for various confidence levels
+_____________________________________________________________
+True Variance = 100.00000
+Difference to detect = 95.16090
+_______________________________________________________________
+Confidence Estimated Estimated
+ Value (%) Sample Size Sample Size
+ (lower one- (upper one-
+ sided test) sided test)
+_______________________________________________________________
+ 50.000 2 2
+ 66.667 2 5
+ 75.000 2 10
+ 90.000 4 32
+ 95.000 5 51
+ 99.000 7 99
+ 99.900 11 174
+ 99.990 15 251
+ 99.999 20 330
+_____________________________________________________________
+Estimated sample sizes required for various confidence levels
+_____________________________________________________________
+True Variance = 100.00000
+Difference to detect = 55.00000
+_______________________________________________________________
+Confidence Estimated Estimated
+ Value (%) Sample Size Sample Size
+ (lower one- (upper one-
+ sided test) sided test)
+_______________________________________________________________
+ 50.000 2 2
+ 66.667 4 10
+ 75.000 8 21
+ 90.000 23 71
+ 95.000 36 115
+ 99.000 71 228
+ 99.900 123 401
+ 99.990 177 580
+ 99.999 232 762
+_____________________________________________________________
+Estimated sample sizes required for various confidence levels
+_____________________________________________________________
+True Variance = 100.00000
+Difference to detect = 1.00000
+_______________________________________________________________
+Confidence Estimated Estimated
+ Value (%) Sample Size Sample Size
+ (lower one- (upper one-
+ sided test) sided test)
+_______________________________________________________________
+ 50.000 2 2
+ 66.667 14696 14993
+ 75.000 36033 36761
+ 90.000 130079 132707
+ 95.000 214283 218612
+ 99.000 428628 437287
+ 99.900 756333 771612
+ 99.990 1095435 1117564
+ 99.999 1440608 1469711
+________________________________________________
+2-Sided Confidence Limits For Standard Deviation
+________________________________________________
+Confidence level (two-sided) = 0.1000000
+Standard Deviation = 1.0000000
+_____________________________________________
+Observations Lower Upper
+ Limit Limit
+_____________________________________________
+ 2 0.5102 15.9472
+ 3 0.5778 4.4154
+ 4 0.6196 2.9200
+ 5 0.6493 2.3724
+ 6 0.6720 2.0893
+ 7 0.6903 1.9154
+ 8 0.7054 1.7972
+ 9 0.7183 1.7110
+ 10 0.7293 1.6452
+ 15 0.7688 1.4597
+ 20 0.7939 1.3704
+ 30 0.8255 1.2797
+ 40 0.8454 1.2320
+ 50 0.8594 1.2017
+ 60 0.8701 1.1805
+ 100 0.8963 1.1336
+ 120 0.9045 1.1203
+ 1000 0.9646 1.0383
+ 10000 0.9885 1.0118
+ 50000 0.9948 1.0052
+ 100000 0.9963 1.0037
+ 1000000 0.9988 1.0012
+________________________________________________
+2-Sided Confidence Limits For Standard Deviation
+________________________________________________
+Confidence level (two-sided) = 0.0500000
+Standard Deviation = 1.0000000
+_____________________________________________
+Observations Lower Upper
+ Limit Limit
+_____________________________________________
+ 2 0.4461 31.9102
+ 3 0.5207 6.2847
+ 4 0.5665 3.7285
+ 5 0.5991 2.8736
+ 6 0.6242 2.4526
+ 7 0.6444 2.2021
+ 8 0.6612 2.0353
+ 9 0.6755 1.9158
+ 10 0.6878 1.8256
+ 15 0.7321 1.5771
+ 20 0.7605 1.4606
+ 30 0.7964 1.3443
+ 40 0.8192 1.2840
+ 50 0.8353 1.2461
+ 60 0.8476 1.2197
+ 100 0.8780 1.1617
+ 120 0.8875 1.1454
+ 1000 0.9580 1.0459
+ 10000 0.9863 1.0141
+ 50000 0.9938 1.0062
+ 100000 0.9956 1.0044
+ 1000000 0.9986 1.0014
+________________________________________________
+2-Sided Confidence Limits For Standard Deviation
+________________________________________________
+Confidence level (two-sided) = 0.0100000
+Standard Deviation = 1.0000000
+_____________________________________________
+Observations Lower Upper
+ Limit Limit
+_____________________________________________
+ 2 0.3562 159.5759
+ 3 0.4344 14.1244
+ 4 0.4834 6.4675
+ 5 0.5188 4.3960
+ 6 0.5464 3.4848
+ 7 0.5688 2.9798
+ 8 0.5875 2.6601
+ 9 0.6036 2.4394
+ 10 0.6177 2.2776
+ 15 0.6686 1.8536
+ 20 0.7018 1.6662
+ 30 0.7444 1.4867
+ 40 0.7718 1.3966
+ 50 0.7914 1.3410
+ 60 0.8065 1.3026
+ 100 0.8440 1.2200
+ 120 0.8558 1.1973
+ 1000 0.9453 1.0609
+ 10000 0.9821 1.0185
+ 50000 0.9919 1.0082
+ 100000 0.9943 1.0058
+ 1000000 0.9982 1.0018
+*/
+
=======================================
--- /dev/null
+++ /trunk/libs/math/example/distribution_construction.cpp Wed Jul 29
07:28:12 2009
@@ -0,0 +1,201 @@
+// distribution_construction.cpp
+
+// Copyright Paul A. Bristow 2007.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Caution: this file contains Quickbook markup as well as code
+// and comments, don't change any of the special comment markups!
+
+//[distribution_construction1
+
+/*`
+
+The structure of distributions is rather different from some other
statistical libraries,
+for example in less object-oriented language like FORTRAN and C,
+that provide a few arguments to each free function.
+This library provides each distribution as a template C++ class.
+A distribution is constructed with a few arguments, and then
+member and non-member functions are used to find values of the
+distribution, often a function of a random variate.
+
+First we need some includes to access the negative binomial distribution
+(and the binomial, beta and gamma too).
+
+*/
+
+#include <boost/math/distributions/negative_binomial.hpp> // for
negative_binomial_distribution
+ using boost::math::negative_binomial_distribution; // default type is
double.
+ using boost::math::negative_binomial; // typedef provides default type
is double.
+#include <boost/math/distributions/binomial.hpp> // for
binomial_distribution.
+#include <boost/math/distributions/beta.hpp> // for beta_distribution.
+#include <boost/math/distributions/gamma.hpp> // for gamma_distribution.
+#include <boost/math/distributions/normal.hpp> // for normal_distribution.
+/*`
+Several examples of constructing distributions follow:
+*/
+//] [/distribution_construction1 end of Quickbook in C++ markup]
+
+
+
+
+int main()
+{
+//[distribution_construction2
+/*`
+First, a negative binomial distribution with 8 successes
+and a success fraction 0.25, 25% or 1 in 4, is constructed like this:
+*/
+ boost::math::negative_binomial_distribution<double> mydist0(8., 0.25);
+ /*`
+ But this is inconveniently long, so by writing
+ */
+ using namespace boost::math;
+ /*`
+ or
+ */
+ using boost::math::negative_binomial_distribution;
+ /*`
+ we can reduce typing.
+
+ Since the vast majority of applications use will be using double
precision,
+ the template argument to the distribution (RealType) defaults
+ to type double, so we can also write:
+ */
+
+ negative_binomial_distribution<> mydist9(8., 0.25); // Uses default
RealType = double.
+
+ /*`
+ But the name "negative_binomial_distribution" is still inconveniently
long,
+ so for most distributions, a convenience typedef is provided, for
example:
+
+ typedef negative_binomial_distribution<double> negative_binomial; //
Reserved name of type double.
+
+ [caution
+ This convenience typedef is /not/ provided if a clash would occur
+ with the name of a function: currently only "beta" and "gamma"
+ fall into this category.
+ ]
+
+ So, after a using statement,
+ */
+
+ using boost::math::negative_binomial;
+
+ /*`
+ we have a convenient typedef to `negative_binomial_distribution<double>`:
+ */
+ negative_binomial mydist(8., 0.25);
+
+ /*`
+ Some more examples using the convenience typedef:
+ */
+ negative_binomial mydist10(5., 0.4); // Both arguments double.
+ /*`
+ And automatic conversion takes place, so you can use integers and floats:
+ */
+ negative_binomial mydist11(5, 0.4); // Using provided typedef double,
int and double arguments.
+ /*`
+ This is probably the most common usage.
+ */
+ negative_binomial mydist12(5., 0.4F); // Double and float arguments.
+ negative_binomial mydist13(5, 1); // Both arguments integer.
+
+ /*`
+ Similarly for most other distributions like the binomial.
+ */
+ binomial mybinomial(1, 0.5); // is more concise than
+ binomial_distribution<> mybinomd1(1, 0.5);
+
+ /*`
+ For cases when the typdef distribution name would clash with a math
special function
+ (currently only beta and gamma)
+ the typedef is deliberately not provided, and the longer version of the
name
+ must be used. For example do not use:
+
+ using boost::math::beta;
+ beta mybetad0(1, 0.5); // Error beta is a math FUNCTION!
+
+ Which produces the error messages:
+
+ [pre
+ error C2146: syntax error : missing ';' before identifier 'mybetad0'
+ warning C4551: function call missing argument list
+ error C3861: 'mybetad0': identifier not found
+ ]
+
+ Instead you should use:
+ */
+ using boost::math::beta_distribution;
+ beta_distribution<> mybetad1(1, 0.5);
+ /*`
+ or for the gamma distribution:
+ */
+ gamma_distribution<> mygammad1(1, 0.5);
+
+ /*`
+ We can, of course, still provide the type explicitly thus:
+ */
+
+ // Explicit double precision:
+ negative_binomial_distribution<double> mydist1(8., 0.25);
+
+ // Explicit float precision, double arguments are truncated to float:
+ negative_binomial_distribution<float> mydist2(8., 0.25);
+
+ // Explicit float precision, integer & double arguments converted to
float.
+ negative_binomial_distribution<float> mydist3(8, 0.25);
+
+ // Explicit float precision, float arguments, so no conversion:
+ negative_binomial_distribution<float> mydist4(8.F, 0.25F);
+
+ // Explicit float precision, integer arguments promoted to float.
+ negative_binomial_distribution<float> mydist5(8, 1);
+
+ // Explicit double precision:
+ negative_binomial_distribution<double> mydist6(8., 0.25);
+
+ // Explicit long double precision:
+ negative_binomial_distribution<long double> mydist7(8., 0.25);
+
+ /*`
+ And if you have your own RealType called MyFPType,
+ for example NTL RR (an arbitrary precision type), then we can write:
+
+ negative_binomial_distribution<MyFPType> mydist6(8, 1); // Integer
arguments -> MyFPType.
+
+ [heading Default arguments to distribution constructors.]
+
+ Note that default constructor arguments are only provided for some
distributions.
+ So if you wrongly assume a default argument you will get an error
message, for example:
+
+ negative_binomial_distribution<> mydist8;
+
+ [pre error C2512 no appropriate default constructor available.]
+
+ No default constructors are provided for the negative binomial,
+ because it is difficult to chose any sensible default values for this
distribution.
+ For other distributions, like the normal distribution,
+ it is obviously very useful to provide 'standard'
+ defaults for the mean and standard deviation thus:
+
+ normal_distribution(RealType mean = 0, RealType sd = 1);
+
+ So in this case we can write:
+ */
+ using boost::math::normal;
+
+ normal norm1; // Standard normal distribution.
+ normal norm2(2); // Mean = 2, std deviation = 1.
+ normal norm3(2, 3); // Mean = 2, std deviation = 3.
+
+ return 0;
+} // int main()
+
+/*`There is no useful output from this program, of course. */
+
+//] [/end of distribution_construction2]
+
=======================================
--- /dev/null
+++ /trunk/libs/math/example/error_handling_example.cpp Wed Jul 29 07:28:12
2009
@@ -0,0 +1,143 @@
+// example_error_handling.cpp
+
+// Copyright Paul A. Bristow 2007.
+// Copyright John Maddock 2007.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Note that this file contains quickbook markup as well as code
+// and comments, don't change any of the special comment markups!
+
+//[error_handling_example
+/*`
+The following example demonstrates the effect of
+setting the macro BOOST_MATH_DOMAIN_ERROR_POLICY
+when an invalid argument is encountered. For the
+purposes of this example, we'll pass a negative
+degrees of freedom parameter to the student's t
+distribution.
+
+Since we know that this is a single file program we could
+just add:
+
+ #define BOOST_MATH_DOMAIN_ERROR_POLICY ignore_error
+
+to the top of the source file to change the default policy
+to one that simply returns a NaN when a domain error occurs.
+Alternatively we could use:
+
+ #define BOOST_MATH_DOMAIN_ERROR_POLICY errno_on_error
+
+To ensure the `::errno` is set when a domain error occurs
+as well as returning a NaN.
+
+This is safe provided the program consists of a single
+translation unit /and/ we place the define /before/ any
+#includes. Note that should we add the define after the includes
+then it will have no effect! A warning such as:
+
+[pre warning C4005: 'BOOST_MATH_OVERFLOW_ERROR_POLICY' : macro
redefinition]
+
+is a certain sign that it will /not/ have the desired effect.
+
+We'll begin our sample program with the needed includes:
+*/
+
+// Boost
+#include <boost/math/distributions/students_t.hpp>
+ using boost::math::students_t; // Probability of students_t(df, t).
+
+// std
+#include <iostream>
+ using std::cout;
+ using std::endl;
+
+#include <stdexcept>
+ using std::exception;
+
+/*`
+
+Next we'll define the program's main() to call the student's t
+distribution with an invalid degrees of freedom parameter, the program
+is set up to handle either an exception or a NaN:
+
+*/
+
+int main()
+{
+ cout << "Example error handling using Student's t function. " << endl;
+ cout << "BOOST_MATH_DOMAIN_ERROR_POLICY is set to: "
+ << BOOST_STRINGIZE(BOOST_MATH_DOMAIN_ERROR_POLICY) << endl;
+
+ double degrees_of_freedom = -1; // A bad argument!
+ double t = 10;
+
+ try
+ {
+ errno = 0;
+ students_t dist(degrees_of_freedom); // exception is thrown here if
enabled
+ double p = cdf(dist, t);
+ // test for error reported by other means:
+ if((boost::math::isnan)(p))
+ {
+ cout << "cdf returned a NaN!" << endl;
+ cout << "errno is set to: " << errno << endl;
+ }
+ else
+ cout << "Probability of Student's t is " << p << endl;
+ }
+ catch(const std::exception& e)
+ {
+ std::cout <<
+ "\n""Message from thrown exception was:\n " << e.what() <<
std::endl;
+ }
+
+ return 0;
+} // int main()
+
+/*`
+
+Here's what the program output looks like with a default build
+(one that does throw exceptions):
+
+[pre
+Example error handling using Student's t function.
+BOOST_MATH_DOMAIN_ERROR_POLICY is set to: throw_on_error
+
+Message from thrown exception was:
+ Error in function
boost::math::students_t_distribution<double>::students_t_distribution:
+ Degrees of freedom argument is -1, but must be > 0 !
+]
+
+Alternatively let's build with:
+
+ #define BOOST_MATH_DOMAIN_ERROR_POLICY ignore_error
+
+Now the program output is:
+
+[pre
+Example error handling using Student's t function.
+BOOST_MATH_DOMAIN_ERROR_POLICY is set to: ignore_error
+cdf returned a NaN!
+errno is set to: 0
+]
+
+And finally let's build with:
+
+ #define BOOST_MATH_DOMAIN_ERROR_POLICY errno_on_error
+
+Which gives the output:
+
+[pre
+Example error handling using Student's t function.
+BOOST_MATH_DOMAIN_ERROR_POLICY is set to: errno_on_error
+cdf returned a NaN!
+errno is set to: 33
+]
+
+*/
+
+//] [error_handling_eg end quickbook markup]
=======================================
--- /dev/null
+++ /trunk/libs/math/example/error_policies_example.cpp Wed Jul 29 07:28:12
2009
@@ -0,0 +1,97 @@
+// error_policies_example.cpp
+
+// Copyright Paul A. Bristow 2007.
+// Copyright John Maddock 2007.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include <boost/math/distributions/normal.hpp>
+ using boost::math::normal_distribution;
+
+#include <boost/math/distributions/students_t.hpp>
+ using boost::math::students_t; // Probability of students_t(df, t).
+ using boost::math::students_t_distribution;
+
+// using namespace boost::math;
+//.\error_policy_normal.cpp(30) : error C2872: 'policy' : ambiguous symbol
+// could
be 'I:\Boost-sandbox\math_toolkit\boost/math/policies/policy.hpp(392) :
boost::math::policies::policy'
+// or 'boost::math::policies'
+
+ // So can't use this using namespace command.
+// Suppose we want a statistical distribution to return infinities,
+// rather than throw exceptions (the default policy), then we can use:
+
+// std
+#include <iostream>
+ using std::cout;
+ using std::endl;
+
+using boost::math::policies::policy;
+// Possible errors
+using boost::math::policies::overflow_error;
+using boost::math::policies::underflow_error;
+using boost::math::policies::domain_error;
+using boost::math::policies::pole_error;
+using boost::math::policies::denorm_error;
+using boost::math::policies::evaluation_error;
+
+using boost::math::policies::ignore_error;
+
+// Define a custom policy to ignore just overflow:
+typedef policy<
+overflow_error<ignore_error>
+ > my_policy;
+
+// Define another custom policy (perhaps ill-advised?)
+// to ignore all errors: domain, pole, overflow, underflow, denorm &
evaluation:
+typedef policy<
+domain_error<ignore_error>,
+pole_error<ignore_error>,
+overflow_error<ignore_error>,
+underflow_error<ignore_error>,
+denorm_error<ignore_error>,
+evaluation_error<ignore_error>
+ > my_ignoreall_policy;
+
+// Define a new distribution with a custom policy to ignore_error
+// (& thus perhaps return infinity for some arguments):
+typedef boost::math::normal_distribution<double, my_policy> my_normal;
+// Note: uses default parameters zero mean and unit standard deviation.
+
+// We could also do the same for another distribution, for example:
+using boost::math::students_t_distribution;
+typedef students_t_distribution<double, my_ignoreall_policy> my_students_t;
+
+int main()
+{
+ cout << "quantile(my_normal(), 0.05); = " << quantile(my_normal(), 0.05)
<< endl; // 0.05 is argument within normal range.
+ cout << "quantile(my_normal(), 0.); = " << quantile(my_normal(), 0.) <<
endl; // argument zero, so expect infinity.
+ cout << "quantile(my_normal(), 0.); = " << quantile(my_normal(), 0.F) <<
endl; // argument zero, so expect infinity.
+
+ cout << "quantile(my_students_t(), 0.); = " <<
quantile(my_students_t(-1), 0.F) << endl; // 'bad' argument negative, so
expect NaN.
+
+ // Construct a (0, 1) normal distribution that ignores all errors,
+ // returning NaN, infinity, zero, or best guess,
+ // and NOT setting errno.
+ normal_distribution<long double, my_ignoreall_policy> my_normal2(0.L,
1.L); // explicit parameters for distribution.
+ cout << "quantile(my_normal2(), 0.); = " << quantile(my_normal2, 0.01)
<< endl; // argument 0.01, so result finite.
+ cout << "quantile(my_normal2(), 0.); = " << quantile(my_normal2, 0.) <<
endl; // argument zero, so expect infinity.
+
+ return 0;
+}
+
+/*
+
+Output:
+
+quantile(my_normal(), 0.05); = -1.64485
+quantile(my_normal(), 0.); = -1.#INF
+quantile(my_normal(), 0.); = -1.#INF
+quantile(my_students_t(), 0.); = 1.#QNAN
+quantile(my_normal2(), 0.); = -2.32635
+quantile(my_normal2(), 0.); = -1.#INF
+
+*/
=======================================
--- /dev/null
+++ /trunk/libs/math/example/error_policy_example.cpp Wed Jul 29 07:28:12
2009
@@ -0,0 +1,88 @@
+// example_policy_handling.cpp
+
+// Copyright Paul A. Bristow 2007.
+// Copyright John Maddock 2007.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// See error_handling_example.cpp for use of
+// macro definition to change policy for
+// domain_error - negative degrees of freedom argument
+// for student's t distribution CDF,
+// and catching the exception.
+
+// See error_handling_policies for more examples.
+
+// Boost
+#include <boost/math/distributions/students_t.hpp>
+using boost::math::students_t_distribution; // Probability of
students_t(df, t).
+using boost::math::students_t; // Probability of students_t(df, t)
convenience typedef for double.
+
+// std
+#include <iostream>
+ using std::cout;
+ using std::endl;
+
+#include <stdexcept>
+ using std::exception;
+
+using boost::math::policies::policy;
+using boost::math::policies::domain_error;
+using boost::math::policies::ignore_error;
+
+// Define a (bad?) policy to ignore domain errors ('bad' arguments):
+typedef policy<
+ domain_error<ignore_error>
+ > my_policy;
+
+// Define my distribution with this different domain error policy:
+typedef students_t_distribution<double, my_policy> my_students_t;
+
+int main()
+{ // Example of error handling of bad argument(s) to a distribution.
+ cout << "Example error handling using Student's t function. " << endl;
+
+ double degrees_of_freedom = -1; double t = -1.; // Two 'bad' arguments!
+
+ try
+ {
+ cout << "Probability of ignore_error Student's t is " <<
cdf(my_students_t(degrees_of_freedom), t) << endl;
+ cout << "Probability of default error policy Student's t is " << endl;
+ // BY contrast the students_t distribution default domain error policy
is to throw,
+ cout << cdf(students_t(-1), -1) << endl; // so this will throw.
+ /*`
+ Message from thrown exception was:
+ Error in function
boost::math::students_t_distribution<double>::students_t_distribution:
+ Degrees of freedom argument is -1, but must be > 0 !
+ */
+
+ // We could also define a 'custom' distribution
+ // with an "ignore overflow error policy" in a single statement:
+ using boost::math::policies::overflow_error;
+ students_t_distribution<double, policy<overflow_error<ignore_error> >
students_t_no_throw(-1);
+
+ }
+ catch(const std::exception& e)
+ {
+ std::cout <<
+ "\n""Message from thrown exception was:\n " << e.what() <<
std::endl;
+ }
+
+ return 0;
+} // int main()
+
+/*
+
+Output:
+
+Example error handling using Student's t function.
+Probability of ignore_error Student's t is 1.#QNAN
+Probability of default error policy Student's t is
+Message from thrown exception was:
+ Error in function
boost::math::students_t_distribution<double>::students_t_distribution:
+ Degrees of freedom argument is -1, but must be > 0 !
+
+*/
=======================================
--- /dev/null
+++ /trunk/libs/math/example/f_test.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,184 @@
+// Copyright John Maddock 2006
+// Copyright Paul A. Bristow 2007
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifdef _MSC_VER
+# pragma warning(disable: 4512) // assignment operator could not be
generated.
+# pragma warning(disable: 4510) // default constructor could not be
generated.
+# pragma warning(disable: 4610) // can never be instantiated - user
defined constructor required.
+#endif
+
+#include <iostream>
+#include <iomanip>
+#include <boost/math/distributions/fisher_f.hpp>
+
+void f_test(
+ double sd1, // Sample 1 std deviation
+ double sd2, // Sample 2 std deviation
+ double N1, // Sample 1 size
+ double N2, // Sample 2 size
+ double alpha) // Significance level
+{
+ //
+ // An F test applied to two sets of data.
+ // We are testing the null hypothesis that the
+ // standard deviation of the samples is equal, and
+ // that any variation is down to chance. We can
+ // also test the alternative hypothesis that any
+ // difference is not down to chance.
+ // See http://www.itl.nist.gov/div898/handbook/eda/section3/eda359.htm
+ //
+ using namespace std;
+ using namespace boost::math;
+
+ // Print header:
+ cout <<
+ "____________________________________\n"
+ "F test for equal standard deviations\n"
+ "____________________________________\n\n";
+ cout << setprecision(5);
+ cout << "Sample 1:\n";
+ cout << setw(55) << left << "Number of Observations" << "= " << N1
<< "\n";
+ cout << setw(55) << left << "Sample Standard Deviation" << "= " << sd1
<< "\n\n";
+ cout << "Sample 2:\n";
+ cout << setw(55) << left << "Number of Observations" << "= " << N2
<< "\n";
+ cout << setw(55) << left << "Sample Standard Deviation" << "= " << sd2
<< "\n\n";
+ //
+ // Now we can calculate and output some stats:
+ //
+ // F-statistic:
+ double F = (sd1 / sd2);
+ F *= F;
+ cout << setw(55) << left << "Test Statistic" << "= " << F << "\n\n";
+ //
+ // Finally define our distribution, and get the probability:
+ //
+ fisher_f dist(N1 - 1, N2 - 1);
+ double p = cdf(dist, F);
+ cout << setw(55) << left << "CDF of test statistic: " << "= "
+ << setprecision(3) << scientific << p << "\n";
+ double ucv = quantile(complement(dist, alpha));
+ double ucv2 = quantile(complement(dist, alpha / 2));
+ double lcv = quantile(dist, alpha);
+ double lcv2 = quantile(dist, alpha / 2);
+ cout << setw(55) << left << "Upper Critical Value at alpha: " << "= "
+ << setprecision(3) << scientific << ucv << "\n";
+ cout << setw(55) << left << "Upper Critical Value at alpha/2: " << "= "
+ << setprecision(3) << scientific << ucv2 << "\n";
+ cout << setw(55) << left << "Lower Critical Value at alpha: " << "= "
+ << setprecision(3) << scientific << lcv << "\n";
+ cout << setw(55) << left << "Lower Critical Value at alpha/2: " << "= "
+ << setprecision(3) << scientific << lcv2 << "\n\n";
+ //
+ // Finally print out results of null and alternative hypothesis:
+ //
+ cout << setw(55) << left <<
+ "Results for Alternative Hypothesis and alpha" << "= "
+ << setprecision(4) << fixed << alpha << "\n\n";
+ cout << "Alternative Hypothesis
Conclusion\n";
+ cout << "Standard deviations are unequal (two sided test) ";
+ if((ucv2 < F) || (lcv2 > F))
+ cout << "NOT REJECTED\n";
+ else
+ cout << "REJECTED\n";
+ cout << "Standard deviation 1 is less than standard deviation 2 ";
+ if(lcv > F)
+ cout << "NOT REJECTED\n";
+ else
+ cout << "REJECTED\n";
+ cout << "Standard deviation 1 is greater than standard deviation 2 ";
+ if(ucv < F)
+ cout << "NOT REJECTED\n";
+ else
+ cout << "REJECTED\n";
+ cout << endl << endl;
+}
+
+int main()
+{
+ //
+ // Run tests for ceramic strength data:
+ // see http://www.itl.nist.gov/div898/handbook/eda/section4/eda42a1.htm
+ // The data for this case study were collected by Said Jahanmir of the
+ // NIST Ceramics Division in 1996 in connection with a NIST/industry
+ // ceramics consortium for strength optimization of ceramic strength.
+ //
+ f_test(65.54909, 61.85425, 240, 240, 0.05);
+ //
+ // And again for the process change comparison:
+ // see http://www.itl.nist.gov/div898/handbook/prc/section3/prc32.htm
+ // A new procedure to assemble a device is introduced and tested for
+ // possible improvement in time of assembly. The question being
addressed
+ // is whether the standard deviation of the new assembly process
(sample 2) is
+ // better (i.e., smaller) than the standard deviation for the old
assembly
+ // process (sample 1).
+ //
+ f_test(4.9082, 2.5874, 11, 9, 0.05);
+ return 0;
+}
+
+/*
+
+Output:
+
+____________________________________
+F test for equal standard deviations
+____________________________________
+
+Sample 1:
+Number of Observations = 240
+Sample Standard Deviation = 65.549
+
+Sample 2:
+Number of Observations = 240
+Sample Standard Deviation = 61.854
+
+Test Statistic = 1.123
+
+CDF of test statistic: = 8.148e-001
+Upper Critical Value at alpha: = 1.238e+000
+Upper Critical Value at alpha/2: = 1.289e+000
+Lower Critical Value at alpha: = 8.080e-001
+Lower Critical Value at alpha/2: = 7.756e-001
+
+Results for Alternative Hypothesis and alpha = 0.0500
+
+Alternative Hypothesis Conclusion
+Standard deviations are unequal (two sided test) REJECTED
+Standard deviation 1 is less than standard deviation 2 REJECTED
+Standard deviation 1 is greater than standard deviation 2 REJECTED
+
+
+____________________________________
+F test for equal standard deviations
+____________________________________
+
+Sample 1:
+Number of Observations = 11.00000
+Sample Standard Deviation = 4.90820
+
+Sample 2:
+Number of Observations = 9.00000
+Sample Standard Deviation = 2.58740
+
+Test Statistic = 3.59847
+
+CDF of test statistic: = 9.589e-001
+Upper Critical Value at alpha: = 3.347e+000
+Upper Critical Value at alpha/2: = 4.295e+000
+Lower Critical Value at alpha: = 3.256e-001
+Lower Critical Value at alpha/2: = 2.594e-001
+
+Results for Alternative Hypothesis and alpha = 0.0500
+
+Alternative Hypothesis Conclusion
+Standard deviations are unequal (two sided test) REJECTED
+Standard deviation 1 is less than standard deviation 2 REJECTED
+Standard deviation 1 is greater than standard deviation 2 NOT REJECTED
+
+*/
+
=======================================
--- /dev/null
+++ /trunk/libs/math/example/find_location_example.cpp Wed Jul 29 07:28:12
2009
@@ -0,0 +1,170 @@
+// find_location.cpp
+
+// Copyright Paul A. Bristow 2007.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Example of finding location (mean)
+// for normal (Gaussian) & Cauchy distribution.
+
+// Note that this file contains Quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+//[find_location1
+/*`
+First we need some includes to access the normal distribution,
+the algorithms to find location (and some std output of course).
+*/
+
+#include <boost/math/distributions/normal.hpp> // for normal_distribution
+ using boost::math::normal; // typedef provides default type is double.
+#include <boost/math/distributions/cauchy.hpp> // for cauchy_distribution
+ using boost::math::cauchy; // typedef provides default type is double.
+#include <boost/math/distributions/find_location.hpp>
+ using boost::math::find_location;
+ using boost::math::complement; // Needed if you want to use the
complement version.
+ using boost::math::policies::policy;
+
+#include <iostream>
+ using std::cout; using std::endl;
+#include <iomanip>
+ using std::setw; using std::setprecision;
+#include <limits>
+ using std::numeric_limits;
+//] [/find_location1]
+
+int main()
+{
+ cout << "Example: Find location (mean)." << endl;
+ try
+ {
+//[find_location2
+/*`
+For this example, we will use the standard normal distribution,
+with mean (location) zero and standard deviation (scale) unity.
+This is also the default for this implementation.
+*/
+ normal N01; // Default 'standard' normal distribution with zero mean and
+ double sd = 1.; // normal default standard deviation is 1.
+/*`Suppose we want to find a different normal distribution whose mean is
shifted
+so that only fraction p (here 0.001 or 0.1%) are below a certain chosen
limit
+(here -2, two standard deviations).
+*/
+ double z = -2.; // z to give prob p
+ double p = 0.001; // only 0.1% below z
+
+ cout << "Normal distribution with mean = " << N01.location()
+ << ", standard deviation " << N01.scale()
+ << ", has " << "fraction <= " << z
+ << ", p = " << cdf(N01, z) << endl;
+ cout << "Normal distribution with mean = " << N01.location()
+ << ", standard deviation " << N01.scale()
+ << ", has " << "fraction > " << z
+ << ", p = " << cdf(complement(N01, z)) << endl; // Note: uses
complement.
+/*`
+[pre
+Normal distribution with mean = 0, standard deviation 1, has fraction <=
-2, p = 0.0227501
+Normal distribution with mean = 0, standard deviation 1, has fraction >
-2, p = 0.97725
+]
+We can now use ''find_location'' to give a new offset mean.
+*/
+ double l = find_location<normal>(z, p, sd);
+ cout << "offset location (mean) = " << l << endl;
+/*`
+that outputs:
+[pre
+offset location (mean) = 1.09023
+]
+showing that we need to shift the mean just over one standard deviation
from its previous value of zero.
+
+Then we can check that we have achieved our objective
+by constructing a new distribution
+with the offset mean (but same standard deviation):
+*/
+ normal np001pc(l, sd); // Same standard_deviation (scale) but with mean
(location) shifted.
+/*`
+And re-calculating the fraction below our chosen limit.
+*/
+cout << "Normal distribution with mean = " << l
+ << " has " << "fraction <= " << z
+ << ", p = " << cdf(np001pc, z) << endl;
+ cout << "Normal distribution with mean = " << l
+ << " has " << "fraction > " << z
+ << ", p = " << cdf(complement(np001pc, z)) << endl;
+/*`
+[pre
+Normal distribution with mean = 1.09023 has fraction <= -2, p = 0.001
+Normal distribution with mean = 1.09023 has fraction > -2, p = 0.999
+]
+
+[h4 Controlling Error Handling from find_location]
+We can also control the policy for handling various errors.
+For example, we can define a new (possibly unwise)
+policy to ignore domain errors ('bad' arguments).
+
+Unless we are using the boost::math namespace, we will need:
+*/
+ using boost::math::policies::policy;
+ using boost::math::policies::domain_error;
+ using boost::math::policies::ignore_error;
+
+/*`
+Using a typedef is often convenient, especially if it is re-used,
+although it is not required, as the various examples below show.
+*/
+ typedef policy<domain_error<ignore_error> > ignore_domain_policy;
+ // find_location with new policy, using typedef.
+ l = find_location<normal>(z, p, sd, ignore_domain_policy());
+ // Default policy policy<>, needs "using boost::math::policies::policy;"
+ l = find_location<normal>(z, p, sd, policy<>());
+ // Default policy, fully specified.
+ l = find_location<normal>(z, p, sd, boost::math::policies::policy<>());
+ // A new policy, ignoring domain errors, without using a typedef.
+ l = find_location<normal>(z, p, sd, policy<domain_error<ignore_error>
());
+/*`
+If we want to use a probability that is the
+[link complements complement of our probability],
+we should not even think of writing `find_location<normal>(z, 1 - p, sd)`,
+but, [link why_complements to avoid loss of accuracy], use the complement
version.
+*/
+ z = 2.;
+ double q = 0.95; // = 1 - p; // complement.
+ l = find_location<normal>(complement(z, q, sd));
+
+ normal np95pc(l, sd); // Same standard_deviation (scale) but with
mean(location) shifted
+ cout << "Normal distribution with mean = " << l << " has "
+ << "fraction <= " << z << " = " << cdf(np95pc, z) << endl;
+ cout << "Normal distribution with mean = " << l << " has "
+ << "fraction > " << z << " = " << cdf(complement(np95pc, z)) << endl;
+ //] [/find_location2]
+ }
+ catch(const std::exception& e)
+ { // Always useful to include try & catch blocks because default policies
+ // are to throw exceptions on arguments that cause errors like
underflow, overflow.
+ // Lacking try & catch blocks, the program will abort without a
message below,
+ // which may give some helpful clues as to the cause of the exception.
+ std::cout <<
+ "\n""Message from thrown exception was:\n " << e.what() <<
std::endl;
+ }
+ return 0;
+} // int main()
+
+
+
+//[find_location_example_output
+/*`
+[pre
+Example: Find location (mean).
+Normal distribution with mean = 0, standard deviation 1, has fraction <=
-2, p = 0.0227501
+Normal distribution with mean = 0, standard deviation 1, has fraction >
-2, p = 0.97725
+offset location (mean) = 1.09023
+Normal distribution with mean = 1.09023 has fraction <= -2, p = 0.001
+Normal distribution with mean = 1.09023 has fraction > -2, p = 0.999
+Normal distribution with mean = 0.355146 has fraction <= 2 = 0.95
+Normal distribution with mean = 0.355146 has fraction > 2 = 0.05
+]
+*/
+//] [/find_location_example_output]
=======================================
--- /dev/null
+++ /trunk/libs/math/example/find_mean_and_sd_normal.cpp Wed Jul 29
07:28:12 2009
@@ -0,0 +1,408 @@
+// find_mean_and_sd_normal.cpp
+
+// Copyright Paul A. Bristow 2007.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Example of finding mean or sd for normal distribution.
+
+// Note that this file contains Quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+//[normal_std
+/*`
+First we need some includes to access the normal distribution,
+the algorithms to find location and scale
+(and some std output of course).
+*/
+
+#include <boost/math/distributions/normal.hpp> // for normal_distribution
+ using boost::math::normal; // typedef provides default type is double.
+#include <boost/math/distributions/cauchy.hpp> // for cauchy_distribution
+ using boost::math::cauchy; // typedef provides default type is double.
+#include <boost/math/distributions/find_location.hpp>
+ using boost::math::find_location;
+#include <boost/math/distributions/find_scale.hpp>
+ using boost::math::find_scale;
+ using boost::math::complement;
+ using boost::math::policies::policy;
+
+#include <iostream>
+ using std::cout; using std::endl; using std::left; using std::showpoint;
using std::noshowpoint;
+#include <iomanip>
+ using std::setw; using std::setprecision;
+#include <limits>
+ using std::numeric_limits;
+//] [/normal_std Quickbook]
+
+int main()
+{
+ cout << "Find_location (mean) and find_scale (standard deviation)
examples." << endl;
+ try
+ {
+
+//[normal_find_location_and_scale_eg
+
+/*`
+[h4 Using find_location and find_scale to meet dispensing and measurement
specifications]
+
+Consider an example from K Krishnamoorthy,
+Handbook of Statistical Distributions with Applications,
+ISBN 1-58488-635-8, (2006) p 126, example 10.3.7.
+
+"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.
+Assume the packing is normally distributed."
+
+We start by constructing a normal distribution with the given parameters:
+*/
+
+double mean = 3.; // kg
+double standard_deviation = 0.1; // kg
+normal packs(mean, standard_deviation);
+/*`We can then find the fraction (or %) of packages that weigh more than
3.1 kg.
+*/
+
+double max_weight = 3.1; // kg
+cout << "Percentage of packs > " << max_weight << " is "
+<< cdf(complement(packs, max_weight)) * 100. << endl; // P(X > 3.1)
+
+/*`We might want to ensure that 95% of packs are over a minimum weight
specification,
+then we want the value of the mean such that P(X < 2.9) = 0.05.
+
+Using the mean of 3 kg, we can estimate
+the fraction of packs that fail to meet the specification of 2.9 kg.
+*/
+
+double minimum_weight = 2.9;
+cout <<"Fraction of packs <= " << minimum_weight << " with a mean of " <<
mean
+ << " is " << cdf(complement(packs, minimum_weight)) << endl;
+// fraction of packs <= 2.9 with a mean of 3 is 0.841345
+
+/*`This is 0.84 - more than the target fraction of 0.95.
+If we want 95% to be over the minimum weight,
+what should we set the mean weight to be?
+
+Using the KK StatCalc program supplied with the book and the method given
on page 126 gives 3.06449.
+
+We can confirm this by constructing a new distribution which we
call 'xpacks'
+with a safety margin mean of 3.06449 thus:
+*/
+double over_mean = 3.06449;
+normal xpacks(over_mean, standard_deviation);
+cout << "Fraction of packs >= " << minimum_weight
+<< " with a mean of " << xpacks.mean()
+ << " is " << cdf(complement(xpacks, minimum_weight)) << endl;
+// fraction of packs >= 2.9 with a mean of 3.06449 is 0.950005
+
+/*`Using this Math Toolkit, we can calculate the required mean directly
thus:
+*/
+double under_fraction = 0.05; // so 95% are above the minimum weight mean
- sd = 2.9
+double low_limit = standard_deviation;
+double offset = mean - low_limit - quantile(packs, under_fraction);
+double nominal_mean = mean + offset;
+// mean + (mean - low_limit - quantile(packs, under_fraction));
+
+normal nominal_packs(nominal_mean, standard_deviation);
+cout << "Setting the packer to " << nominal_mean << " will mean that "
+ << "fraction of packs >= " << minimum_weight
+ << " is " << cdf(complement(nominal_packs, minimum_weight)) << endl;
+// Setting the packer to 3.06449 will mean that fraction of packs >= 2.9
is 0.95
+
+/*`
+This calculation is generalized as the free function called
+[link math_toolkit.dist.dist_ref.dist_algorithms find_location].
+
+To use this we will need to
+*/
+
+#include <boost/math/distributions/find_location.hpp>
+ using boost::math::find_location;
+/*`and then use find_location function to find safe_mean,
+& construct a new normal distribution called 'goodpacks'.
+*/
+double safe_mean = find_location<normal>(minimum_weight, under_fraction,
standard_deviation);
+normal good_packs(safe_mean, standard_deviation);
+/*`with the same confirmation as before:
+*/
+cout << "Setting the packer to " << nominal_mean << " will mean that "
+ << "fraction of packs >= " << minimum_weight
+ << " is " << cdf(complement(good_packs, minimum_weight)) << endl;
+// Setting the packer to 3.06449 will mean that fraction of packs >= 2.9
is 0.95
+
+/*`
+[h4 Using Cauchy-Lorentz instead of normal distribution]
+
+After examining the weight distribution of a large number of packs, we
might decide that,
+after all, the assumption of a normal distribution is not really justified.
+We might find that the fit is better to a __cauchy_distrib.
+This distribution has wider 'wings', so that whereas most of the values
+are closer to the mean than the normal, there are also more values
than 'normal'
+that lie further from the mean than the normal.
+
+This might happen because a larger than normal lump of meat is either
included or excluded.
+
+We first create a __cauchy_distrib with the original mean and standard
deviation,
+and estimate the fraction that lie below our minimum weight specification.
+*/
+
+cauchy cpacks(mean, standard_deviation);
+cout << "Cauchy Setting the packer to " << mean << " will mean that "
+ << "fraction of packs >= " << minimum_weight
+ << " is " << cdf(complement(cpacks, minimum_weight)) << endl;
+// Cauchy Setting the packer to 3 will mean that fraction of packs >= 2.9
is 0.75
+
+/*`Note that far fewer of the packs meet the specification, only 75%
instead of 95%.
+Now we can repeat the find_location, using the cauchy distribution as
template parameter,
+in place of the normal used above.
+*/
+
+double lc = find_location<cauchy>(minimum_weight, under_fraction,
standard_deviation);
+cout << "find_location<cauchy>(minimum_weight, over fraction,
standard_deviation); " << lc << endl;
+// find_location<cauchy>(minimum_weight, over fraction,
packs.standard_deviation()); 3.53138
+/*`Note that the safe_mean setting needs to be much higher, 3.53138
instead of 3.06449,
+so we will make rather less profit.
+
+And again confirm that the fraction meeting specification is as expected.
+*/
+cauchy goodcpacks(lc, standard_deviation);
+cout << "Cauchy Setting the packer to " << lc << " will mean that "
+ << "fraction of packs >= " << minimum_weight
+ << " is " << cdf(complement(goodcpacks, minimum_weight)) << endl;
+// Cauchy Setting the packer to 3.53138 will mean that fraction of packs
= 2.9 is 0.95
+
+/*`Finally we could estimate the effect of a much tighter specification,
+that 99% of packs met the specification.
+*/
+
+cout << "Cauchy Setting the packer to "
+ << find_location<cauchy>(minimum_weight, 0.99, standard_deviation)
+ << " will mean that "
+ << "fraction of packs >= " << minimum_weight
+ << " is " << cdf(complement(goodcpacks, minimum_weight)) << endl;
+
+/*`Setting the packer to 3.13263 will mean that fraction of packs >= 2.9
is 0.99,
+but will more than double the mean loss from 0.0644 to 0.133 kg per pack.
+
+Of course, this calculation is not limited to packs of meat, it applies to
dispensing anything,
+and it also applies to a 'virtual' material like any measurement.
+
+The only caveat is that the calculation assumes that the standard
deviation (scale) is known with
+a reasonably low uncertainty, something that is not so easy to ensure in
practice.
+And that the distribution is well defined, __normal_distrib or
__cauchy_distrib, or some other.
+
+If one is simply dispensing a very large number of packs,
+then it may be feasible to measure the weight of hundreds or thousands of
packs.
+With a healthy 'degrees of freedom', the confidence intervals for the
standard deviation
+are not too wide, typically about + and - 10% for hundreds of observations.
+
+For other applications, where it is more difficult or expensive to make
many observations,
+the confidence intervals are depressingly wide.
+
+See [link math_toolkit.dist.stat_tut.weg.cs_eg.chi_sq_intervals Confidence
Intervals on the standard deviation]
+for a worked example
+[@../../../example/chi_square_std_dev_test.cpp chi_square_std_dev_test.cpp]
+of estimating these intervals.
+
+
+[h4 Changing the scale or standard deviation]
+
+Alternatively, we could invest in a better (more precise) packer
+(or measuring device) with a lower standard deviation, or scale.
+
+This might cost more, but would reduce the amount we have to 'give away'
+in order to meet the specification.
+
+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
+*/
+double p = 0.05; // wanted p th quantile.
+cout << "Quantile of " << p << " = " << quantile(packs, p)
+ << ", mean = " << packs.mean() << ", sd = " <<
packs.standard_deviation() << endl;
+/*`
+Quantile of 0.05 = 2.83551, mean = 3, sd = 0.1
+
+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.
+
+Let's start by guessing that it (now 0.1) needs to be halved, to a
standard deviation of 0.05 kg.
+*/
+normal pack05(mean, 0.05);
+cout << "Quantile of " << p << " = " << quantile(pack05, p)
+ << ", mean = " << pack05.mean() << ", sd = " <<
pack05.standard_deviation() << endl;
+// Quantile of 0.05 = 2.91776, mean = 3, sd = 0.05
+
+cout <<"Fraction of packs >= " << minimum_weight << " with a mean of " <<
mean
+ << " and standard deviation of " << pack05.standard_deviation()
+ << " is " << cdf(complement(pack05, minimum_weight)) << endl;
+// Fraction of packs >= 2.9 with a mean of 3 and standard deviation of
0.05 is 0.97725
+/*`
+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 standard deviation to 0.06
kg,
+constructing another new distribution called pack06.
+*/
+normal pack06(mean, 0.06);
+cout << "Quantile of " << p << " = " << quantile(pack06, p)
+ << ", mean = " << pack06.mean() << ", sd = " <<
pack06.standard_deviation() << endl;
+// Quantile of 0.05 = 2.90131, mean = 3, sd = 0.06
+
+cout <<"Fraction of packs >= " << minimum_weight << " with a mean of " <<
mean
+ << " and standard deviation of " << pack06.standard_deviation()
+ << " is " << cdf(complement(pack06, minimum_weight)) << endl;
+// Fraction of packs >= 2.9 with a mean of 3 and standard deviation of
0.06 is 0.95221
+/*`
+Now we are getting really close, but to do the job properly,
+we might need to use root finding method, for example the tools provided,
+and used elsewhere, in the Math Toolkit, see
+[link math_toolkit.toolkit.internals1.roots2 Root Finding Without
Derivatives].
+
+But in this (normal) distribution case, we can and should be even smarter
+and make a direct calculation.
+*/
+
+/*`Our required limit is minimum_weight = 2.9 kg, often called the random
variate z.
+For a standard normal distribution, then probability p = N((minimum_weight
- mean) / sd).
+
+We want to find the standard deviation that would be required to meet this
limit,
+so that the p th quantile is located at z (minimum_weight).
+In this case, the 0.05 (5%) quantile is at 2.9 kg pack weight, when the
mean is 3 kg,
+ensuring that 0.95 (95%) of packs are above the minimum weight.
+
+Rearranging, we can directly calculate the required standard deviation:
+*/
+normal N01; // standard normal distribution with meamn zero and unit
standard deviation.
+p = 0.05;
+double qp = quantile(N01, p);
+double sd95 = (minimum_weight - mean) / qp;
+
+cout << "For the "<< p << "th quantile to be located at "
+ << minimum_weight << ", would need a standard deviation of " << sd95 <<
endl;
+// For the 0.05th quantile to be located at 2.9, would need a standard
deviation of 0.0607957
+
+/*`We can now construct a new (normal) distribution pack95 for
the 'better' packer,
+and check that our distribution will meet the specification.
+*/
+
+normal pack95(mean, sd95);
+cout <<"Fraction of packs >= " << minimum_weight << " with a mean of " <<
mean
+ << " and standard deviation of " << pack95.standard_deviation()
+ << " is " << cdf(complement(pack95, minimum_weight)) << endl;
+// Fraction of packs >= 2.9 with a mean of 3 and standard deviation of
0.0607957 is 0.95
+
+/*`This calculation is generalized in the free function find_scale,
+as shown below, giving the same standard deviation.
+*/
+double ss = find_scale<normal>(minimum_weight, under_fraction,
packs.mean());
+cout << "find_scale<normal>(minimum_weight, under_fraction,
packs.mean()); " << ss << endl;
+// find_scale<normal>(minimum_weight, under_fraction, packs.mean());
0.0607957
+
+/*`If we had defined an over_fraction, or percentage that must pass
specification
+*/
+double over_fraction = 0.95;
+/*`And (wrongly) written
+
+ double sso = find_scale<normal>(minimum_weight, over_fraction,
packs.mean());
+
+With the default policy, we would get a message like
+
+[pre
+Message from thrown exception was:
+ Error in function boost::math::find_scale<Dist, Policy>(double, double,
double, Policy):
+ Computed scale (-0.060795683191176959) is <= 0! Was the complement
intended?
+]
+
+But this would return a *negative* standard deviation - obviously
impossible.
+The probability should be 1 - over_fraction, not over_fraction, thus:
+*/
+
+double ss1o = find_scale<normal>(minimum_weight, 1 - over_fraction,
packs.mean());
+cout << "find_scale<normal>(minimum_weight, under_fraction,
packs.mean()); " << ss1o << endl;
+// find_scale<normal>(minimum_weight, under_fraction, packs.mean());
0.0607957
+
+/*`But notice that using '1 - over_fraction' - will lead to a
+[link why_complements loss of accuracy, especially if over_fraction was
close to unity.]
+In this (very common) case, we should instead use the __complements,
+giving the most accurate result.
+*/
+
+double ssc = find_scale<normal>(complement(minimum_weight, over_fraction,
packs.mean()));
+cout << "find_scale<normal>(complement(minimum_weight, over_fraction,
packs.mean())); " << ssc << endl;
+// find_scale<normal>(complement(minimum_weight, over_fraction,
packs.mean())); 0.0607957
+
+/*`Note that our guess of 0.06 was close to the accurate value of
0.060795683191176959.
+
+We can again confirm our prediction thus:
+*/
+
+normal pack95c(mean, ssc);
+cout <<"Fraction of packs >= " << minimum_weight << " with a mean of " <<
mean
+ << " and standard deviation of " << pack95c.standard_deviation()
+ << " is " << cdf(complement(pack95c, minimum_weight)) << endl;
+// Fraction of packs >= 2.9 with a mean of 3 and standard deviation of
0.0607957 is 0.95
+
+/*`Notice that these two deceptively simple questions:
+
+* Do we over-fill to make sure we meet a minimum specification (or
under-fill to avoid an overdose)?
+
+and/or
+
+* Do we measure better?
+
+are actually extremely common.
+
+The weight of beef might be replaced by a measurement of more or less
anything,
+from drug tablet content, Apollo landing rocket firing, X-ray treatment
doses...
+
+The scale can be variation in dispensing or uncertainty in measurement.
+*/
+//] [/normal_find_location_and_scale_eg Quickbook end]
+
+ }
+ catch(const std::exception& e)
+ { // Always useful to include try & catch blocks because default policies
+ // are to throw exceptions on arguments that cause errors like
underflow, overflow.
+ // Lacking try & catch blocks, the program will abort without a
message below,
+ // which may give some helpful clues as to the cause of the exception.
+ std::cout <<
+ "\n""Message from thrown exception was:\n " << e.what() <<
std::endl;
+ }
+ return 0;
+} // int main()
+
+
+/*
+
+Output is:
+
+Find_location and find_scale examples.
+Percentage of packs > 3.1 is 15.8655
+Fraction of packs <= 2.9 with a mean of 3 is 0.841345
+Fraction of packs >= 2.9 with a mean of 3.06449 is 0.950005
+Setting the packer to 3.06449 will mean that fraction of packs >= 2.9 is
0.95
+Setting the packer to 3.06449 will mean that fraction of packs >= 2.9 is
0.95
+Cauchy Setting the packer to 3 will mean that fraction of packs >= 2.9 is
0.75
+find_location<cauchy>(minimum_weight, over fraction, standard_deviation);
3.53138
+Cauchy Setting the packer to 3.53138 will mean that fraction of packs >=
2.9 is 0.95
+Cauchy Setting the packer to -0.282052 will mean that fraction of packs >=
2.9 is 0.95
+Quantile of 0.05 = 2.83551, mean = 3, sd = 0.1
+Quantile of 0.05 = 2.91776, mean = 3, sd = 0.05
+Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.05
is 0.97725
+Quantile of 0.05 = 2.90131, mean = 3, sd = 0.06
+Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.06
is 0.95221
+For the 0.05th quantile to be located at 2.9, would need a standard
deviation of 0.0607957
+Fraction of packs >= 2.9 with a mean of 3 and standard deviation of
0.0607957 is 0.95
+find_scale<normal>(minimum_weight, under_fraction, packs.mean()); 0.0607957
+find_scale<normal>(minimum_weight, under_fraction, packs.mean()); 0.0607957
+find_scale<normal>(complement(minimum_weight, over_fraction,
packs.mean())); 0.0607957
+Fraction of packs >= 2.9 with a mean of 3 and standard deviation of
0.0607957 is 0.95
+
+*/
+
+
+
=======================================
--- /dev/null
+++ /trunk/libs/math/example/find_root_example.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,239 @@
+// find_root_example.cpp
+
+// Copyright Paul A. Bristow 2007.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Example of using root finding.
+
+// Note that this file contains Quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+//[root_find1
+/*`
+First we need some includes to access the normal distribution
+(and some std output of course).
+*/
+
+#include <boost/math/tools/roots.hpp> // root finding.
+
+#include <boost/math/distributions/normal.hpp> // for normal_distribution
+ using boost::math::normal; // typedef provides default type is double.
+
+#include <iostream>
+ using std::cout; using std::endl; using std::left; using std::showpoint;
using std::noshowpoint;
+#include <iomanip>
+ using std::setw; using std::setprecision;
+#include <limits>
+ using std::numeric_limits;
+//] //[/root_find1]
+
+namespace boost{ namespace math { namespace tools
+{
+
+template <class F, class T, class Tol>
+inline std::pair<T, T> bracket_and_solve_root(F f, // functor
+ const T& guess,
+ const T& factor,
+ bool rising,
+ Tol tol, // binary functor
specifying termination when tol(min, max) becomes true.
+ // eps_tolerance most
suitable for this continuous function
+ boost::uintmax_t& max_iter);
// explicit (rather than default) max iterations.
+// return interval as a pair containing result.
+
+namespace detail
+{
+
+ // Functor for finding standard deviation:
+ template <class RealType, class Policy>
+ struct standard_deviation_functor
+ {
+ standard_deviation_functor(RealType m, RealType s, RealType d)
+ : mean(m), standard_deviation(s)
+ {
+ }
+ RealType operator()(const RealType& sd)
+ { // Unary functor - the function whose root is to be found.
+ if(sd <= tools::min_value<RealType>())
+ { //
+ return 1;
+ }
+ normal_distribution<RealType, Policy> t(mean, sd);
+ RealType qa = quantile(complement(t, alpha));
+ RealType qb = quantile(complement(t, beta));
+ qa += qb;
+ qa *= qa;
+ qa *= ratio;
+ qa -= (df + 1);
+ return qa;
+ } // operator()
+ RealType mean;
+ RealType standard_deviation;
+ }; // struct standard_deviation_functor
+} // namespace detail
+
+template <class RealType, class Policy>
+RealType normal_distribution<RealType, Policy>::find_standard_deviation(
+ RealType difference_from_mean,
+ RealType mean,
+ RealType sd,
+ RealType hint) // Best guess available - current sd if none better?
+{
+ static const char* function
= "boost::math::normal_distribution<%1%>::find_standard_deviation";
+
+ // Check for domain errors:
+ RealType error_result;
+ if(false == detail::check_probability(
+ function, sd, &error_result, Policy())
+ )
+ return error_result;
+
+ if(hint <= 0)
+ { // standard deviation can never be negative.
+ hint = 1;
+ }
+
+ detail::standard_deviation_functor<RealType, Policy> f(mean, sd,
difference_from_mean);
+ tools::eps_tolerance<RealType> tol(policies::digits<RealType,
Policy>());
+ boost::uintmax_t max_iter = 100;
+ std::pair<RealType, RealType> r = tools::bracket_and_solve_root(f,
hint, RealType(2), false, tol, max_iter, Policy());
+ RealType result = r.first + (r.second - r.first) / 2;
+ if(max_iter == 100)
+ {
+ policies::raise_evaluation_error<RealType>(function, "Unable to
locate solution in a reasonable time:"
+ " either there is no answer to how many degrees of freedom are
required"
+ " or the answer is infinite. Current best guess is %1%", result,
Policy());
+ }
+ return result;
+} // find_standard_deviation
+} // namespace tools
+} // namespace math
+} // namespace boost
+
+
+int main()
+{
+ cout << "Example: Normal distribution, root finding.";
+ try
+ {
+
+//[root_find2
+
+/*`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.
+*/
+
+double mean = 3.; // kg
+double standard_deviation = 0.1; // kg
+normal packs(mean, standard_deviation);
+
+double max_weight = 3.1; // kg
+cout << "Percentage of packs > " << max_weight << " is "
+<< cdf(complement(packs, max_weight)) << endl; // P(X > 3.1)
+
+double under_weight = 2.9;
+cout <<"fraction of packs <= " << under_weight << " with a mean of " <<
mean
+ << " is " << cdf(complement(packs, under_weight)) << endl;
+// fraction of packs <= 2.9 with a mean of 3 is 0.841345
+// This is 0.84 - more than the target 0.95
+// Want 95% to be over this weight, so what should we set the mean weight
to be?
+// KK StatCalc says:
+double over_mean = 3.0664;
+normal xpacks(over_mean, standard_deviation);
+cout << "fraction of packs >= " << under_weight
+<< " with a mean of " << xpacks.mean()
+ << " is " << cdf(complement(xpacks, under_weight)) << endl;
+// fraction of packs >= 2.9 with a mean of 3.06449 is 0.950005
+double under_fraction = 0.05; // so 95% are above the minimum weight mean
- sd = 2.9
+double low_limit = standard_deviation;
+double offset = mean - low_limit - quantile(packs, under_fraction);
+double nominal_mean = mean + offset;
+
+normal nominal_packs(nominal_mean, standard_deviation);
+cout << "Setting the packer to " << nominal_mean << " will mean that "
+ << "fraction of packs >= " << under_weight
+ << " is " << cdf(complement(nominal_packs, under_weight)) << endl;
+
+/*`
+Setting the packer to 3.06449 will mean that fraction of packs >= 2.9 is
0.95.
+
+Setting the packer to 3.13263 will mean that fraction of packs >= 2.9 is
0.99,
+but will more than double the mean loss from 0.0644 to 0.133.
+
+Alternatively, we could invest in a better (more precise) packer with a
lower standard deviation.
+
+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
+*/
+double p = 0.05; // wanted p th quantile.
+cout << "Quantile of " << p << " = " << quantile(packs, p)
+ << ", mean = " << packs.mean() << ", sd = " <<
packs.standard_deviation() << endl; //
+/*`
+Quantile of 0.05 = 2.83551, mean = 3, sd = 0.1
+
+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.
+
+Let's start by guessing that it (now 0.1) needs to be halved, to a
standard deviation of 0.05
+*/
+normal pack05(mean, 0.05);
+cout << "Quantile of " << p << " = " << quantile(pack05, p)
+ << ", mean = " << pack05.mean() << ", sd = " <<
pack05.standard_deviation() << endl;
+
+cout <<"Fraction of packs >= " << under_weight << " with a mean of " <<
mean
+ << " and standard deviation of " << pack05.standard_deviation()
+ << " is " << cdf(complement(pack05, under_weight)) << endl;
+//
+/*`
+Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.05
is 0.9772
+
+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
+*/
+
+normal pack06(mean, 0.06);
+cout << "Quantile of " << p << " = " << quantile(pack06, p)
+ << ", mean = " << pack06.mean() << ", sd = " <<
pack06.standard_deviation() << endl;
+
+cout <<"Fraction of packs >= " << under_weight << " with a mean of " <<
mean
+ << " and standard deviation of " << pack06.standard_deviation()
+ << " is " << cdf(complement(pack06, under_weight)) << endl;
+/*`
+Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.06
is 0.9522
+
+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
+[link math_toolkit.toolkit.internals1.roots2 Root Finding Without
Derivatives].
+
+But in this normal distribution case, we could be even smarter and make a
direct calculation.
+*/
+//] [/root_find2]
+
+ }
+ catch(const std::exception& e)
+ { // Always useful to include try & catch blocks because default policies
+ // are to throw exceptions on arguments that cause errors like
underflow, overflow.
+ // Lacking try & catch blocks, the program will abort without a
message below,
+ // which may give some helpful clues as to the cause of the exception.
+ std::cout <<
+ "\n""Message from thrown exception was:\n " << e.what() <<
std::endl;
+ }
+ return 0;
+} // int main()
+
+/*
+
+Output is:
+
+
+
+*/
=======================================
--- /dev/null
+++ /trunk/libs/math/example/find_scale_example.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,180 @@
+// find_scale.cpp
+
+// Copyright Paul A. Bristow 2007.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Example of finding scale (standard deviation) for normal (Gaussian).
+
+// Note that this file contains Quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+//[find_scale1
+/*`
+First we need some includes to access the __normal_distrib,
+the algorithms to find scale (and some std output of course).
+*/
+
+#include <boost/math/distributions/normal.hpp> // for normal_distribution
+ using boost::math::normal; // typedef provides default type is double.
+#include <boost/math/distributions/find_scale.hpp>
+ using boost::math::find_scale;
+ using boost::math::complement; // Needed if you want to use the
complement version.
+ using boost::math::policies::policy; // Needed to specify the error
handling policy.
+
+#include <iostream>
+ using std::cout; using std::endl;
+#include <iomanip>
+ using std::setw; using std::setprecision;
+#include <limits>
+ using std::numeric_limits;
+//] [/find_scale1]
+
+int main()
+{
+ cout << "Example: Find scale (standard deviation)." << endl;
+ try
+ {
+//[find_scale2
+/*`
+For this example, we will use the standard __normal_distrib,
+with location (mean) zero and standard deviation (scale) unity.
+Conveniently, this is also the default for this implementation's
constructor.
+*/
+ normal N01; // Default 'standard' normal distribution with zero mean
+ double sd = 1.; // and standard deviation is 1.
+/*`Suppose we want to find a different normal distribution with standard
deviation
+so that only fraction p (here 0.001 or 0.1%) are below a certain chosen
limit
+(here -2. standard deviations).
+*/
+ double z = -2.; // z to give prob p
+ double p = 0.001; // only 0.1% below z = -2
+
+ cout << "Normal distribution with mean = " << N01.location() // aka
N01.mean()
+ << ", standard deviation " << N01.scale() // aka
N01.standard_deviation()
+ << ", has " << "fraction <= " << z
+ << ", p = " << cdf(N01, z) << endl;
+ cout << "Normal distribution with mean = " << N01.location()
+ << ", standard deviation " << N01.scale()
+ << ", has " << "fraction > " << z
+ << ", p = " << cdf(complement(N01, z)) << endl; // Note: uses
complement.
+/*`
+[pre
+Normal distribution with mean = 0 has fraction <= -2, p = 0.0227501
+Normal distribution with mean = 0 has fraction > -2, p = 0.97725
+]
+Noting that p = 0.02 instead of our target of 0.001,
+we can now use `find_scale` to give a new standard deviation.
+*/
+ double l = N01.location();
+ double s = find_scale<normal>(z, p, l);
+ cout << "scale (standard deviation) = " << s << endl;
+/*`
+that outputs:
+[pre
+scale (standard deviation) = 0.647201
+]
+showing that we need to reduce the standard deviation from 1. to 0.65.
+
+Then we can check that we have achieved our objective
+by constructing a new distribution
+with the new standard deviation (but same zero mean):
+*/
+ normal np001pc(N01.location(), s);
+/*`
+And re-calculating the fraction below (and above) our chosen limit.
+*/
+ cout << "Normal distribution with mean = " << l
+ << " has " << "fraction <= " << z
+ << ", p = " << cdf(np001pc, z) << endl;
+ cout << "Normal distribution with mean = " << l
+ << " has " << "fraction > " << z
+ << ", p = " << cdf(complement(np001pc, z)) << endl;
+/*`
+[pre
+Normal distribution with mean = 0 has fraction <= -2, p = 0.001
+Normal distribution with mean = 0 has fraction > -2, p = 0.999
+]
+
+[h4 Controlling how Errors from find_scale are handled]
+We can also control the policy for handling various errors.
+For example, we can define a new (possibly unwise)
+policy to ignore domain errors ('bad' arguments).
+
+Unless we are using the boost::math namespace, we will need:
+*/
+ using boost::math::policies::policy;
+ using boost::math::policies::domain_error;
+ using boost::math::policies::ignore_error;
+
+/*`
+Using a typedef is convenient, especially if it is re-used,
+although it is not required, as the various examples below show.
+*/
+ typedef policy<domain_error<ignore_error> > ignore_domain_policy;
+ // find_scale with new policy, using typedef.
+ l = find_scale<normal>(z, p, l, ignore_domain_policy());
+ // Default policy policy<>, needs using boost::math::policies::policy;
+
+ l = find_scale<normal>(z, p, l, policy<>());
+ // Default policy, fully specified.
+ l = find_scale<normal>(z, p, l, boost::math::policies::policy<>());
+ // New policy, without typedef.
+ l = find_scale<normal>(z, p, l, policy<domain_error<ignore_error> >());
+/*`
+If we want to express a probability, say 0.999, that is a complement, `1 -
p`
+we should not even think of writing `find_scale<normal>(z, 1 - p, l)`,
+but [link why_complements instead], use the __complements version.
+*/
+ z = -2.;
+ double q = 0.999; // = 1 - p; // complement of 0.001.
+ sd = find_scale<normal>(complement(z, q, l));
+
+ normal np95pc(l, sd); // Same standard_deviation (scale) but with
mean(scale) shifted
+ cout << "Normal distribution with mean = " << l << " has "
+ << "fraction <= " << z << " = " << cdf(np95pc, z) << endl;
+ cout << "Normal distribution with mean = " << l << " has "
+ << "fraction > " << z << " = " << cdf(complement(np95pc, z)) << endl;
+
+/*`
+Sadly, it is all too easy to get probabilities the wrong way round,
+when you may get a warning like this:
+[pre
+Message from thrown exception was:
+ Error in function boost::math::find_scale<Dist,
Policy>(complement(double, double, double, Policy)):
+ Computed scale (-0.48043523852179076) is <= 0! Was the complement
intended?
+]
+The default error handling policy is to throw an exception with this
message,
+but if you chose a policy to ignore the error,
+the (impossible) negative scale is quietly returned.
+*/
+//] [/find_scale2]
+ }
+ catch(const std::exception& e)
+ { // Always useful to include try & catch blocks because default policies
+ // are to throw exceptions on arguments that cause errors like
underflow, overflow.
+ // Lacking try & catch blocks, the program will abort without a
message below,
+ // which may give some helpful clues as to the cause of the exception.
+ std::cout <<
+ "\n""Message from thrown exception was:\n " << e.what() <<
std::endl;
+ }
+ return 0;
+} // int main()
+
+//[find_scale_example_output
+/*`
+[pre
+Example: Find scale (standard deviation).
+Normal distribution with mean = 0, standard deviation 1, has fraction <=
-2, p = 0.0227501
+Normal distribution with mean = 0, standard deviation 1, has fraction >
-2, p = 0.97725
+scale (standard deviation) = 0.647201
+Normal distribution with mean = 0 has fraction <= -2, p = 0.001
+Normal distribution with mean = 0 has fraction > -2, p = 0.999
+Normal distribution with mean = 0.946339 has fraction <= -2 = 0.001
+Normal distribution with mean = 0.946339 has fraction > -2 = 0.999
+]
+*/
+//] [/find_scale_example_output]
=======================================
--- /dev/null
+++ /trunk/libs/math/example/nc_chi_sq_example.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,109 @@
+// Copyright John Maddock 2008
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Caution: this file contains Quickbook markup as well as code
+// and comments, don't change any of the special comment markups!
+
+//[nccs_eg
+
+/*`
+
+This example computes a table of the power of the [chi][super 2]
+test at the 5% significance level, for various degrees of freedom
+and non-centrality parameters. The table is deliberately the same
+as Table 6 from "The Non-Central [chi][super 2] and F-Distributions and
+their applications.", P. B. Patnaik, Biometrika, Vol. 36, No. 1/2 (June
1949),
+202-232.
+
+First we need some includes to access the non-central chi squared
distribution
+(and some basic std output of course).
+
+*/
+
+#include <boost/math/distributions/non_central_chi_squared.hpp>
+#include <iostream>
+
+int main()
+{
+ /*`
+ Create a table of the power of the [chi][super 2] test at
+ 5% significance level, start with a table header:
+ */
+ std::cout << "[table\n[[[nu]]";
+ for(int lam = 2; lam <= 20; lam += 2)
+ {
+ std::cout << "[[lambda]=" << lam << "]";
+ }
+ std::cout << "]\n";
+
+ /*`
+ (Note: the enclosing [] brackets are to format as a table in
Boost.Quickbook).
+
+ Enumerate the rows and columns and print the power of the test
+ for each table cell:
+ */
+
+ for(int n = 2; n <= 20; ++n)
+ {
+ std::cout << "[[" << n << "]";
+ for(int lam = 2; lam <= 20; lam += 2)
+ {
+ /*`
+ Calculate the [chi][super 2] statistic for a 5% significance:
+ */
+ double cs = quantile(complement(boost::math::chi_squared(n),
0.05));
+ /*`
+ The power of the test is given by the complement of the CDF
+ of the non-central [chi][super 2] distribution:
+ */
+ double beta =
cdf(complement(boost::math::non_central_chi_squared(n, lam), cs));
+ /*`
+ Then output the cell value:
+ */
+ std::cout << "[" << std::setprecision(3) << beta << "]";
+ }
+ std::cout << "]" << std::endl;
+ }
+ std::cout << "]" << std::endl;
+}
+
+/*`
+The output from this program is a table in Boost.Quickbook format as shown
below.
+
+We can interpret this as follows - for example if [nu]=10 and [lambda]=10
+then the power of the test is 0.542 - so we have only a 54% chance of
+correctly detecting that our null hypothesis is false, and a 46% chance
+of incurring a type II error (failing to reject the null hypothesis when
+it is in fact false):
+
+[table
+[[[nu]][[lambda]=2][[lambda]=4][[lambda]=6][[lambda]=8][[lambda]=10][[lambda]=12][[lambda]=14][[lambda]=16][[lambda]=18][[lambda]=20]]
+[[2][0.226][0.415][0.584][0.718][0.815][0.883][0.928][0.957][0.974][0.985]]
+[[3][0.192][0.359][0.518][0.654][0.761][0.84][0.896][0.934][0.959][0.975]]
+[[4][0.171][0.32][0.47][0.605][0.716][0.802][0.866][0.912][0.943][0.964]]
+[[5][0.157][0.292][0.433][0.564][0.677][0.769][0.839][0.89][0.927][0.952]]
+[[6][0.146][0.27][0.403][0.531][0.644][0.738][0.813][0.869][0.911][0.94]]
+[[7][0.138][0.252][0.378][0.502][0.614][0.71][0.788][0.849][0.895][0.928]]
+[[8][0.131][0.238][0.357][0.477][0.588][0.685][0.765][0.829][0.879][0.915]]
+[[9][0.125][0.225][0.339][0.454][0.564][0.661][0.744][0.811][0.863][0.903]]
+[[10][0.121][0.215][0.323][0.435][0.542][0.64][0.723][0.793][0.848][0.891]]
+[[11][0.117][0.206][0.309][0.417][0.523][0.62][0.704][0.775][0.833][0.878]]
+[[12][0.113][0.198][0.297][0.402][0.505][0.601][0.686][0.759][0.818][0.866]]
+[[13][0.11][0.191][0.286][0.387][0.488][0.584][0.669][0.743][0.804][0.854]]
+[[14][0.108][0.185][0.276][0.374][0.473][0.567][0.653][0.728][0.791][0.842]]
+[[15][0.105][0.179][0.267][0.362][0.459][0.552][0.638][0.713][0.777][0.83]]
+[[16][0.103][0.174][0.259][0.351][0.446][0.538][0.623][0.699][0.764][0.819]]
+[[17][0.101][0.169][0.251][0.341][0.434][0.525][0.609][0.686][0.752][0.807]]
+[[18][0.0992][0.165][0.244][0.332][0.423][0.512][0.596][0.673][0.74][0.796]]
+[[19][0.0976][0.161][0.238][0.323][0.412][0.5][0.584][0.66][0.728][0.786]]
+[[20][0.0961][0.158][0.232][0.315][0.402][0.489][0.572][0.648][0.716][0.775]]
+]
+
+See [@../../../example/nc_chi_sq_example.cpp nc_chi_sq_example.cpp] for
the full C++ source code.
+
+*/
+
+//]
=======================================
--- /dev/null
+++ /trunk/libs/math/example/neg_binom_confidence_limits.cpp Wed Jul 29
07:28:12 2009
@@ -0,0 +1,179 @@
+// neg_binomial_confidence_limits.cpp
+
+// Copyright John Maddock 2006
+// Copyright Paul A. Bristow 2007
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Caution: this file contains quickbook markup as well as code
+// and comments, don't change any of the special comment markups!
+
+//[neg_binomial_confidence_limits
+
+/*`
+
+First we need some includes to access the negative binomial distribution
+(and some basic std output of course).
+
+*/
+
+#include <boost/math/distributions/negative_binomial.hpp>
+using boost::math::negative_binomial;
+
+#include <iostream>
+using std::cout; using std::endl;
+#include <iomanip>
+using std::setprecision;
+using std::setw; using std::left; using std::fixed; using std::right;
+
+/*`
+First define a table of significance levels: these are the
+probabilities that the true occurrence frequency lies outside the
calculated
+interval:
+*/
+
+ double alpha[] = { 0.5, 0.25, 0.1, 0.05, 0.01, 0.001, 0.0001, 0.00001 };
+
+/*`
+
+Confidence value as % is (1 - alpha) * 100, so alpha 0.05 == 95% confidence
+that the true occurence frequency lies *inside* the calculated interval.
+
+We need a function to calculate and print confidence limits
+for an observed frequency of occurrence
+that follows a negative binomial distribution.
+
+*/
+
+void confidence_limits_on_frequency(unsigned trials, unsigned successes)
+{
+ // trials = Total number of trials.
+ // successes = Total number of observed successes.
+ // failures = trials - successes.
+ // success_fraction = successes /trials.
+ // Print out general info:
+ cout <<
+ "______________________________________________\n"
+ "2-Sided Confidence Limits For Success Fraction\n"
+ "______________________________________________\n\n";
+ cout << setprecision(7);
+ cout << setw(40) << left << "Number of trials" << " = " << trials
<< "\n";
+ cout << setw(40) << left << "Number of successes" << " = " <<
successes << "\n";
+ cout << setw(40) << left << "Number of failures" << " = " << trials -
successes << "\n";
+ cout << setw(40) << left << "Observed frequency of occurrence" << "
= " << double(successes) / trials << "\n";
+
+ // Print table header:
+ cout << "\n\n"
+ "___________________________________________\n"
+ "Confidence Lower Upper\n"
+ " Value (%) Limit Limit\n"
+ "___________________________________________\n";
+
+
+/*`
+And now for the important part - the bounds themselves.
+For each value of /alpha/, we call `find_lower_bound_on_p` and
+`find_upper_bound_on_p` to obtain lower and upper bounds respectively.
+Note that since we are calculating a two-sided interval,
+we must divide the value of alpha in two. Had we been calculating a
+single-sided interval, for example: ['"Calculate a lower bound so that we
are P%
+sure that the true occurrence frequency is greater than some value"]
+then we would *not* have divided by two.
+
+*/
+
+ // Now print out the upper and lower limits for the alpha table values.
+ for(unsigned i = 0; i < sizeof(alpha)/sizeof(alpha[0]); ++i)
+ {
+ // Confidence value:
+ cout << fixed << setprecision(3) << setw(10) << right << 100 *
(1-alpha[i]);
+ // Calculate bounds:
+ double lower = negative_binomial::find_lower_bound_on_p(trials,
successes, alpha[i]/2);
+ double upper = negative_binomial::find_upper_bound_on_p(trials,
successes, alpha[i]/2);
+ // Print limits:
+ cout << fixed << setprecision(5) << setw(15) << right << lower;
+ cout << fixed << setprecision(5) << setw(15) << right << upper <<
endl;
+ }
+ cout << endl;
+} // void confidence_limits_on_frequency(unsigned trials, unsigned
successes)
+
+/*`
+
+And then call confidence_limits_on_frequency with increasing numbers of
trials,
+but always the same success fraction 0.1, or 1 in 10.
+
+*/
+
+int main()
+{
+ confidence_limits_on_frequency(20, 2); // 20 trials, 2 successes, 2 in
20, = 1 in 10 = 0.1 success fraction.
+ confidence_limits_on_frequency(200, 20); // More trials, but same 0.1
success fraction.
+ confidence_limits_on_frequency(2000, 200); // Many more trials, but same
0.1 success fraction.
+
+ return 0;
+} // int main()
+
+//] [/negative_binomial_confidence_limits_eg end of Quickbook in C++
markup]
+
+/*
+
+______________________________________________
+2-Sided Confidence Limits For Success Fraction
+______________________________________________
+Number of trials = 20
+Number of successes = 2
+Number of failures = 18
+Observed frequency of occurrence = 0.1
+___________________________________________
+Confidence Lower Upper
+ Value (%) Limit Limit
+___________________________________________
+ 50.000 0.04812 0.13554
+ 75.000 0.03078 0.17727
+ 90.000 0.01807 0.22637
+ 95.000 0.01235 0.26028
+ 99.000 0.00530 0.33111
+ 99.900 0.00164 0.41802
+ 99.990 0.00051 0.49202
+ 99.999 0.00016 0.55574
+______________________________________________
+2-Sided Confidence Limits For Success Fraction
+______________________________________________
+Number of trials = 200
+Number of successes = 20
+Number of failures = 180
+Observed frequency of occurrence = 0.1000000
+___________________________________________
+Confidence Lower Upper
+ Value (%) Limit Limit
+___________________________________________
+ 50.000 0.08462 0.11350
+ 75.000 0.07580 0.12469
+ 90.000 0.06726 0.13695
+ 95.000 0.06216 0.14508
+ 99.000 0.05293 0.16170
+ 99.900 0.04343 0.18212
+ 99.990 0.03641 0.20017
+ 99.999 0.03095 0.21664
+______________________________________________
+2-Sided Confidence Limits For Success Fraction
+______________________________________________
+Number of trials = 2000
+Number of successes = 200
+Number of failures = 1800
+Observed frequency of occurrence = 0.1000000
+___________________________________________
+Confidence Lower Upper
+ Value (%) Limit Limit
+___________________________________________
+ 50.000 0.09536 0.10445
+ 75.000 0.09228 0.10776
+ 90.000 0.08916 0.11125
+ 95.000 0.08720 0.11352
+ 99.000 0.08344 0.11802
+ 99.900 0.07921 0.12336
+ 99.990 0.07577 0.12795
+ 99.999 0.07282 0.13206
+*/
=======================================
--- /dev/null
+++ /trunk/libs/math/example/neg_binomial_sample_sizes.cpp Wed Jul 29
07:28:12 2009
@@ -0,0 +1,202 @@
+// neg_binomial_sample_sizes.cpp
+
+// Copyright Paul A. Bristow 2007
+// Copyright John Maddock 2006
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include <boost/math/distributions/negative_binomial.hpp>
+using boost::math::negative_binomial;
+
+// Default RealType is double so this permits use of:
+double find_minimum_number_of_trials(
+double k, // number of failures (events), k >= 0.
+double p, // fraction of trails for which event occurs, 0 <= p <= 1.
+double probability); // probability threshold, 0 <= probability <= 1.
+
+#include <iostream>
+using std::cout;
+using std::endl;
+using std::fixed;
+using std::right;
+#include <iomanip>
+using std::setprecision;
+using std::setw;
+
+//[neg_binomial_sample_sizes
+
+/*`
+It centres around a routine that prints out
+a table of minimum sample sizes for various probability thresholds:
+*/
+ void find_number_of_trials(double failures, double p);
+/*`
+First define a table of significance levels: these are the maximum
+acceptable probability that /failure/ or fewer events will be observed.
+*/
+ double alpha[] = { 0.5, 0.25, 0.1, 0.05, 0.01, 0.001, 0.0001, 0.00001 };
+/*`
+Confidence value as % is (1 - alpha) * 100, so alpha 0.05 == 95% confidence
+that the desired number of failures will be observed.
+
+Much of the rest of the program is pretty-printing, the important part
+is in the calculation of minimum number of trials required for each
+value of alpha using:
+
+ (int)ceil(negative_binomial::find_minimum_number_of_trials(failures, p,
alpha[i]);
+
+*/
+
+/*`
+
+find_minimum_number_of_trials returns a double,
+so ceil rounds this up to ensure we have an integral minimum number of
trials.
+
+*/
+
+void find_number_of_trials(double failures, double p)
+{
+ // trials = number of trials
+ // failures = number of failures before achieving required success(es).
+ // p = success fraction (0 <= p <= 1.).
+ //
+ // Calculate how many trials we need to ensure the
+ // required number of failures DOES exceed "failures".
+
+ cout << "\n""Target number of failures = " << failures;
+ cout << ", Success fraction = " << 100 * p << "%" << endl;
+ // Print table header:
+ cout << "\n\n"
+ "____________________________\n"
+ "Confidence Min Number\n"
+ " Value (%) Of Trials \n"
+ "____________________________\n";
+ // Now print out the data for the alpha table values.
+ for(unsigned i = 0; i < sizeof(alpha)/sizeof(alpha[0]); ++i)
+ { // Confidence values %:
+ cout << fixed << setprecision(3) << setw(10) << right << 100 *
(1-alpha[i]) << " "
+ // find_minimum_number_of_trials
+ << setw(6) << right
+ <<
(int)ceil(negative_binomial::find_minimum_number_of_trials(failures, p,
alpha[i]))
+ << endl;
+ }
+ cout << endl;
+} // void find_number_of_trials(double failures, double p)
+
+/*` finally we can produce some tables of minimum trials for the chosen
confidence levels:
+*/
+
+int main()
+{
+ find_number_of_trials(5, 0.5);
+ find_number_of_trials(50, 0.5);
+ find_number_of_trials(500, 0.5);
+ find_number_of_trials(50, 0.1);
+ find_number_of_trials(500, 0.1);
+ find_number_of_trials(5, 0.9);
+
+ return 0;
+} // int main()
+
+//] [/neg_binomial_sample_sizes.cpp end of Quickbook in C++ markup]
+
+/*
+
+Output is:
+
+Autorun
"i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\neg_binomial_sample_sizes.exe"
+Target number of failures = 5, Success fraction = 50%
+____________________________
+Confidence Min Number
+ Value (%) Of Trials
+____________________________
+ 50.000 11
+ 75.000 14
+ 90.000 17
+ 95.000 18
+ 99.000 22
+ 99.900 27
+ 99.990 31
+ 99.999 36
+Target number of failures = 50.000, Success fraction = 50.000%
+____________________________
+Confidence Min Number
+ Value (%) Of Trials
+____________________________
+ 50.000 101
+ 75.000 109
+ 90.000 115
+ 95.000 119
+ 99.000 128
+ 99.900 137
+ 99.990 146
+ 99.999 154
+Target number of failures = 500.000, Success fraction = 50.000%
+____________________________
+Confidence Min Number
+ Value (%) Of Trials
+____________________________
+ 50.000 1001
+ 75.000 1023
+ 90.000 1043
+ 95.000 1055
+ 99.000 1078
+ 99.900 1104
+ 99.990 1126
+ 99.999 1146
+Target number of failures = 50.000, Success fraction = 10.000%
+____________________________
+Confidence Min Number
+ Value (%) Of Trials
+____________________________
+ 50.000 56
+ 75.000 58
+ 90.000 60
+ 95.000 61
+ 99.000 63
+ 99.900 66
+ 99.990 68
+ 99.999 71
+Target number of failures = 500.000, Success fraction = 10.000%
+____________________________
+Confidence Min Number
+ Value (%) Of Trials
+____________________________
+ 50.000 556
+ 75.000 562
+ 90.000 567
+ 95.000 570
+ 99.000 576
+ 99.900 583
+ 99.990 588
+ 99.999 594
+Target number of failures = 5.000, Success fraction = 90.000%
+____________________________
+Confidence Min Number
+ Value (%) Of Trials
+____________________________
+ 50.000 57
+ 75.000 73
+ 90.000 91
+ 95.000 103
+ 99.000 127
+ 99.900 159
+ 99.990 189
+ 99.999 217
+Target number of failures = 5.000, Success fraction = 40.000%
+____________________________
+Confidence Min Number
+ Value (%) Of Trials
+____________________________
+ 50.000 10
+ 75.000 11
+ 90.000 13
+ 95.000 15
+ 99.000 18
+ 99.900 21
+ 99.990 25
+ 99.999 28
+
+*/
=======================================
--- /dev/null
+++ /trunk/libs/math/example/negative_binomial_example1.cpp Wed Jul 29
07:28:12 2009
@@ -0,0 +1,521 @@
+// negative_binomial_example1.cpp
+
+// Copyright Paul A. Bristow 2007.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Example 1 of using negative_binomial distribution.
+
+//[negative_binomial_eg1_1
+
+/*`
+Based on [@http://en.wikipedia.org/wiki/Negative_binomial_distribution
+a problem by Dr. Diane Evans,
+Professor of Mathematics at Rose-Hulman Institute of Technology].
+
+Pat is required to sell candy bars to raise money for the 6th grade field
trip.
+There are thirty houses in the neighborhood,
+and Pat is not supposed to return home until five candy bars have been
sold.
+So the child goes door to door, selling candy bars.
+At each house, there is a 0.4 probability (40%) of selling one candy bar
+and a 0.6 probability (60%) of selling nothing.
+
+What is the probability mass (density) function (pdf) for selling the last
(fifth)
+candy bar at the nth house?
+
+The Negative Binomial(r, p) distribution describes the probability of k
failures
+and r successes in k+r Bernoulli(p) trials with success on the last trial.
+(A [@http://en.wikipedia.org/wiki/Bernoulli_distribution Bernoulli trial]
+is one with only two possible outcomes, success of failure,
+and p is the probability of success).
+See also [@ http://en.wikipedia.org/wiki/Bernoulli_distribution Bernoulli
distribution]
+and [@http://www.math.uah.edu/stat/bernoulli/Introduction.xhtml Bernoulli
applications].
+
+In this example, we will deliberately produce a variety of calculations
+and outputs to demonstrate the ways that the negative binomial distribution
+can be implemented with this library: it is also deliberately
over-commented.
+
+First we need to #define macros to control the error and discrete handling
policies.
+For this simple example, we want to avoid throwing
+an exception (the default policy) and just return infinity.
+We want to treat the distribution as if it was continuous,
+so we choose a discrete_quantile policy of real,
+rather than the default policy integer_round_outwards.
+*/
+#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error
+#define BOOST_MATH_DISCRETE_QUANTILE_POLICY real
+/*`
+After that we need some includes to provide easy access to the negative
binomial distribution,
+[caution It is vital to #include distributions etc *after* the above
#defines]
+and we need some std library iostream, of course.
+*/
+#include <boost/math/distributions/negative_binomial.hpp>
+ // for negative_binomial_distribution
+ using boost::math::negative_binomial; // typedef provides default type
is double.
+ using ::boost::math::pdf; // Probability mass function.
+ using ::boost::math::cdf; // Cumulative density function.
+ using ::boost::math::quantile;
+
+#include <iostream>
+ using std::cout; using std::endl;
+ using std::noshowpoint; using std::fixed; using std::right; using
std::left;
+#include <iomanip>
+ using std::setprecision; using std::setw;
+
+#include <limits>
+ using std::numeric_limits;
+//] [negative_binomial_eg1_1]
+
+int main()
+{
+ cout <<"Selling candy bars - using the negative binomial distribution."
+ << "\nby Dr. Diane Evans,"
+ "\nProfessor of Mathematics at Rose-Hulman Institute of Technology,"
+ << "\nsee
http://en.wikipedia.org/wiki/Negative_binomial_distribution\n";
+ << endl;
+ cout << endl;
+ cout.precision(5);
+ // None of the values calculated have a useful accuracy as great this,
but
+ // INF shows wrongly with < 5 !
+ //
https://connect.microsoft.com/VisualStudio/feedback/ViewFeedback.aspx?FeedbackID=240227
+//[negative_binomial_eg1_2
+/*`
+It is always sensible to use try and catch blocks because defaults
policies are to
+throw an exception if anything goes wrong.
+
+A simple catch block (see below) will ensure that you get a
+helpful error message instead of an abrupt program abort.
+*/
+ try
+ {
+/*`
+Selling five candy bars means getting five successes, so successes r = 5.
+The total number of trials (n, in this case, houses visited) this takes is
therefore
+ = sucesses + failures or k + r = k + 5.
+*/
+ double sales_quota = 5; // Pat's sales quota - successes (r).
+/*`
+At each house, there is a 0.4 probability (40%) of selling one candy bar
+and a 0.6 probability (60%) of selling nothing.
+*/
+ double success_fraction = 0.4; // success_fraction (p) - so
failure_fraction is 0.6.
+/*`
+The Negative Binomial(r, p) distribution describes the probability of k
failures
+and r successes in k+r Bernoulli(p) trials with success on the last trial.
+(A [@http://en.wikipedia.org/wiki/Bernoulli_distribution Bernoulli trial]
+is one with only two possible outcomes, success of failure,
+and p is the probability of success).
+
+We therefore start by constructing a negative binomial distribution
+with parameters sales_quota (required successes) and probability of
success.
+*/
+ negative_binomial nb(sales_quota, success_fraction); // type double by
default.
+/*`
+To confirm, display the success_fraction & successes parameters of the
distribution.
+*/
+ cout << "Pat has a sales per house success rate of " <<
success_fraction
+ << ".\nTherefore he would, on average, sell " <<
nb.success_fraction() * 100
+ << " bars after trying 100 houses." << endl;
+
+ int all_houses = 30; // The number of houses on the estate.
+
+ cout << "With a success rate of " << nb.success_fraction()
+ << ", he might expect, on average,\n"
+ "to need to visit about " << success_fraction * all_houses
+ << " houses in order to sell all " << nb.successes() << " bars. "
<< endl;
+/*`
+[pre
+Pat has a sales per house success rate of 0.4.
+Therefore he would, on average, sell 40 bars after trying 100 houses.
+With a success rate of 0.4, he might expect, on average,
+to need to visit about 12 houses in order to sell all 5 bars.
+]
+
+The random variable of interest is the number of houses
+that must be visited to sell five candy bars,
+so we substitute k = n - 5 into a negative_binomial(5, 0.4)
+and obtain the [link math.dist.pdf probability mass (density) function
(pdf or pmf)]
+of the distribution of houses visited.
+Obviously, the best possible case is that Pat makes sales on all the first
five houses.
+
+We calculate this using the pdf function:
+*/
+ cout << "Probability that Pat finishes on the " << sales_quota << "th
house is "
+ << pdf(nb, 5 - sales_quota) << endl; // == pdf(nb, 0)
+/*`
+Of course, he could not finish on fewer than 5 houses because he must sell
5 candy bars.
+So the 5th house is the first that he could possibly finish on.
+
+To finish on or before the 8th house, Pat must finish at the 5th, 6th, 7th
or 8th house.
+The probability that he will finish on *exactly* ( == ) on any house
+is the Probability Density Function (pdf).
+*/
+ cout << "Probability that Pat finishes on the 6th house is "
+ << pdf(nb, 6 - sales_quota) << endl;
+ cout << "Probability that Pat finishes on the 7th house is "
+ << pdf(nb, 7 - sales_quota) << endl;
+ cout << "Probability that Pat finishes on the 8th house is "
+ << pdf(nb, 8 - sales_quota) << endl;
+/*`
+[pre
+Probability that Pat finishes on the 6th house is 0.03072
+Probability that Pat finishes on the 7th house is 0.055296
+Probability that Pat finishes on the 8th house is 0.077414
+]
+
+The sum of the probabilities for these houses is the Cumulative
Distribution Function (cdf).
+We can calculate it by adding the individual probabilities.
+*/
+ cout << "Probability that Pat finishes on or before the 8th house is
sum "
+ "\n" << "pdf(sales_quota) + pdf(6) + pdf(7) + pdf(8) = "
+ // Sum each of the mass/density probabilities for houses sales_quota
= 5, 6, 7, & 8.
+ << pdf(nb, 5 - sales_quota) // 0 failures.
+ + pdf(nb, 6 - sales_quota) // 1 failure.
+ + pdf(nb, 7 - sales_quota) // 2 failures.
+ + pdf(nb, 8 - sales_quota) // 3 failures.
+ << endl;
+/*`[pre
+pdf(sales_quota) + pdf(6) + pdf(7) + pdf(8) = 0.17367
+]
+
+Or, usually better, by using the negative binomial *cumulative*
distribution function.
+*/
+ cout << "\nProbability of selling his quota of " << sales_quota
+ << " bars\non or before the " << 8 << "th house is "
+ << cdf(nb, 8 - sales_quota) << endl;
+/*`[pre
+Probability of selling his quota of 5 bars on or before the 8th house is
0.17367
+]*/
+ cout << "\nProbability that Pat finishes exactly on the 10th house is "
+ << pdf(nb, 10 - sales_quota) << endl;
+ cout << "\nProbability of selling his quota of " << sales_quota
+ << " bars\non or before the " << 10 << "th house is "
+ << cdf(nb, 10 - sales_quota) << endl;
+/*`
+[pre
+Probability that Pat finishes exactly on the 10th house is 0.10033
+Probability of selling his quota of 5 bars on or before the 10th house is
0.3669
+]*/
+ cout << "Probability that Pat finishes exactly on the 11th house is "
+ << pdf(nb, 11 - sales_quota) << endl;
+ cout << "\nProbability of selling his quota of " << sales_quota
+ << " bars\non or before the " << 11 << "th house is "
+ << cdf(nb, 11 - sales_quota) << endl;
+/*`[pre
+Probability that Pat finishes on the 11th house is 0.10033
+Probability of selling his quota of 5 candy bars
+on or before the 11th house is 0.46723
+]*/
+ cout << "Probability that Pat finishes exactly on the 12th house is "
+ << pdf(nb, 12 - sales_quota) << endl;
+
+ cout << "\nProbability of selling his quota of " << sales_quota
+ << " bars\non or before the " << 12 << "th house is "
+ << cdf(nb, 12 - sales_quota) << endl;
+/*`[pre
+Probability that Pat finishes on the 12th house is 0.094596
+Probability of selling his quota of 5 candy bars
+on or before the 12th house is 0.56182
+]
+Finally consider the risk of Pat not selling his quota of 5 bars
+even after visiting all the houses.
+Calculate the probability that he /will/ sell on
+or before the last house:
+Calculate the probability that he would sell all his quota on the very
last house.
+*/
+ cout << "Probability that Pat finishes on the " << all_houses
+ << " house is " << pdf(nb, all_houses - sales_quota) << endl;
+/*`
+Probability of selling his quota of 5 bars on the 30th house is
+[pre
+Probability that Pat finishes on the 30 house is 0.00069145
+]
+when he'd be very unlucky indeed!
+
+What is the probability that Pat exhausts all 30 houses in the
neighborhood,
+and *still* doesn't sell the required 5 candy bars?
+*/
+ cout << "\nProbability of selling his quota of " << sales_quota
+ << " bars\non or before the " << all_houses << "th house is "
+ << cdf(nb, all_houses - sales_quota) << endl;
+/*`
+[pre
+Probability of selling his quota of 5 bars
+on or before the 30th house is 0.99849
+]
+
+/*`So the risk of failing even after visiting all the houses is 1 - this
probability,
+ ``1 - cdf(nb, all_houses - sales_quota``
+But using this expression may cause serious inaccuracy,
+so it would be much better to use the complement of the cdf:
+So the risk of failing even at, or after, the 31th (non-existent) houses
is 1 - this probability,
+ ``1 - cdf(nb, all_houses - sales_quota)``
+But using this expression may cause serious inaccuracy.
+So it would be much better to use the complement of the cdf.
+[link why_complements Why complements?]
+*/
+ cout << "\nProbability of failing to sell his quota of " << sales_quota
+ << " bars\neven after visiting all " << all_houses << " houses is "
+ << cdf(complement(nb, all_houses - sales_quota)) << endl;
+/*`
+[pre
+Probability of failing to sell his quota of 5 bars
+even after visiting all 30 houses is 0.0015101
+]
+We can also use the quantile (percentile), the inverse of the cdf, to
+predict which house Pat will finish on. So for the 8th house:
+*/
+ double p = cdf(nb, (8 - sales_quota));
+ cout << "Probability of meeting sales quota on or before 8th house is "<<
p << endl;
+/*`
+[pre
+Probability of meeting sales quota on or before 8th house is 0.174
+]
+*/
+ cout << "If the confidence of meeting sales quota is " << p
+ << ", then the finishing house is " << quantile(nb, p) +
sales_quota << endl;
+
+ cout<< " quantile(nb, p) = " << quantile(nb, p) << endl;
+/*`
+[pre
+If the confidence of meeting sales quota is 0.17367, then the finishing
house is 8
+]
+Demanding absolute certainty that all 5 will be sold,
+implies an infinite number of trials.
+(Of course, there are only 30 houses on the estate,
+so he can't ever be *certain* of selling his quota).
+*/
+ cout << "If the confidence of meeting sales quota is " << 1.
+ << ", then the finishing house is " << quantile(nb, 1) +
sales_quota << endl;
+ // 1.#INF == infinity.
+/*`[pre
+If the confidence of meeting sales quota is 1, then the finishing house is
1.#INF
+]
+And similarly for a few other probabilities:
+*/
+ cout << "If the confidence of meeting sales quota is " << 0.
+ << ", then the finishing house is " << quantile(nb, 0.) +
sales_quota << endl;
+
+ cout << "If the confidence of meeting sales quota is " << 0.5
+ << ", then the finishing house is " << quantile(nb, 0.5) +
sales_quota << endl;
+
+ cout << "If the confidence of meeting sales quota is " << 1 - 0.00151
// 30 th
+ << ", then the finishing house is " << quantile(nb, 1 - 0.00151) +
sales_quota << endl;
+/*`
+[pre
+If the confidence of meeting sales quota is 0, then the finishing house is
5
+If the confidence of meeting sales quota is 0.5, then the finishing house
is 11.337
+If the confidence of meeting sales quota is 0.99849, then the finishing
house is 30
+]
+
+Notice that because we chose a discrete quantile policy of real,
+the result can be an 'unreal' fractional house.
+
+If the opposite is true, we don't want to assume any confidence, then this
is tantamount
+to assuming that all the first sales_quota trials will be successful sales.
+*/
+ cout << "If confidence of meeting quota is zero\n(we assume all houses
are successful sales)"
+ ", then finishing house is " << sales_quota << endl;
+/*`
+[pre
+If confidence of meeting quota is zero (we assume all houses are
successful sales), then finishing house is 5
+If confidence of meeting quota is 0, then finishing house is 5
+]
+We can list quantiles for a few probabilities:
+*/
+
+ double ps[] = {0., 0.001, 0.01, 0.05, 0.1, 0.5, 0.9, 0.95, 0.99,
0.999, 1.};
+ // Confidence as fraction = 1-alpha, as percent = 100 * (1-alpha[i]) %
+ cout.precision(3);
+ for (int i = 0; i < sizeof(ps)/sizeof(ps[0]); i++)
+ {
+ cout << "If confidence of meeting quota is " << ps[i]
+ << ", then finishing house is " << quantile(nb, ps[i]) +
sales_quota
+ << endl;
+ }
+
+/*`
+[pre
+If confidence of meeting quota is 0, then finishing house is 5
+If confidence of meeting quota is 0.001, then finishing house is 5
+If confidence of meeting quota is 0.01, then finishing house is 5
+If confidence of meeting quota is 0.05, then finishing house is 6.2
+If confidence of meeting quota is 0.1, then finishing house is 7.06
+If confidence of meeting quota is 0.5, then finishing house is 11.3
+If confidence of meeting quota is 0.9, then finishing house is 17.8
+If confidence of meeting quota is 0.95, then finishing house is 20.1
+If confidence of meeting quota is 0.99, then finishing house is 24.8
+If confidence of meeting quota is 0.999, then finishing house is 31.1
+If confidence of meeting quota is 1, then finishing house is 1.#INF
+]
+
+We could have applied a ceil function to obtain a 'worst case' integer
value for house.
+``ceil(quantile(nb, ps[i]))``
+
+Or, if we had used the default discrete quantile policy, integer_outside,
by omitting
+``#define BOOST_MATH_DISCRETE_QUANTILE_POLICY real``
+we would have achieved the same effect.
+
+The real result gives some suggestion which house is most likely.
+For example, compare the real and integer_outside for 95% confidence.
+
+[pre
+If confidence of meeting quota is 0.95, then finishing house is 20.1
+If confidence of meeting quota is 0.95, then finishing house is 21
+]
+The real value 20.1 is much closer to 20 than 21, so integer_outside is
pessimistic.
+We could also use integer_round_nearest policy to suggest that 20 is more
likely.
+
+Finally, we can tabulate the probability for the last sale being exactly
on each house.
+*/
+ cout << "\nHouse for " << sales_quota << "th (last) sale. Probability
(%)" << endl;
+ cout.precision(5);
+ for (int i = (int)sales_quota; i < all_houses+1; i++)
+ {
+ cout << left << setw(3) << i << " " <<
setw(8) << cdf(nb, i - sales_quota) << endl;
+ }
+ cout << endl;
+/*`
+[pre
+House for 5 th (last) sale. Probability (%)
+5 0.01024
+6 0.04096
+7 0.096256
+8 0.17367
+9 0.26657
+10 0.3669
+11 0.46723
+12 0.56182
+13 0.64696
+14 0.72074
+15 0.78272
+16 0.83343
+17 0.874
+18 0.90583
+19 0.93039
+20 0.94905
+21 0.96304
+22 0.97342
+23 0.98103
+24 0.98655
+25 0.99053
+26 0.99337
+27 0.99539
+28 0.99681
+29 0.9978
+30 0.99849
+]
+
+As noted above, using a catch block is always a good idea, even if you do
not expect to use it.
+*/
+ }
+ catch(const std::exception& e)
+ { // Since we have set an overflow policy of ignore_error,
+ // an overflow exception should never be thrown.
+ std::cout << "\nMessage from thrown exception was:\n " << e.what() <<
std::endl;
+/*`
+For example, without a ignore domain error policy, if we asked for
``pdf(nb, -1)`` for example, we would get:
+[pre
+Message from thrown exception was:
+ Error in function boost::math::pdf(const
negative_binomial_distribution<double>&, double):
+ Number of failures argument is -1, but must be >= 0 !
+]
+*/
+//] [/ negative_binomial_eg1_2]
+ }
+ return 0;
+} // int main()
+
+
+/*
+
+Output is:
+
+Selling candy bars - using the negative binomial distribution.
+by Dr. Diane Evans,
+Professor of Mathematics at Rose-Hulman Institute of Technology,
+see http://en.wikipedia.org/wiki/Negative_binomial_distribution
+Pat has a sales per house success rate of 0.4.
+Therefore he would, on average, sell 40 bars after trying 100 houses.
+With a success rate of 0.4, he might expect, on average,
+to need to visit about 12 houses in order to sell all 5 bars.
+Probability that Pat finishes on the 5th house is 0.01024
+Probability that Pat finishes on the 6th house is 0.03072
+Probability that Pat finishes on the 7th house is 0.055296
+Probability that Pat finishes on the 8th house is 0.077414
+Probability that Pat finishes on or before the 8th house is sum
+pdf(sales_quota) + pdf(6) + pdf(7) + pdf(8) = 0.17367
+Probability of selling his quota of 5 bars
+on or before the 8th house is 0.17367
+Probability that Pat finishes exactly on the 10th house is 0.10033
+Probability of selling his quota of 5 bars
+on or before the 10th house is 0.3669
+Probability that Pat finishes exactly on the 11th house is 0.10033
+Probability of selling his quota of 5 bars
+on or before the 11th house is 0.46723
+Probability that Pat finishes exactly on the 12th house is 0.094596
+Probability of selling his quota of 5 bars
+on or before the 12th house is 0.56182
+Probability that Pat finishes on the 30 house is 0.00069145
+Probability of selling his quota of 5 bars
+on or before the 30th house is 0.99849
+Probability of failing to sell his quota of 5 bars
+even after visiting all 30 houses is 0.0015101
+Probability of meeting sales quota on or before 8th house is 0.17367
+If the confidence of meeting sales quota is 0.17367, then the finishing
house is 8
+ quantile(nb, p) = 3
+If the confidence of meeting sales quota is 1, then the finishing house is
1.#INF
+If the confidence of meeting sales quota is 0, then the finishing house is
5
+If the confidence of meeting sales quota is 0.5, then the finishing house
is 11.337
+If the confidence of meeting sales quota is 0.99849, then the finishing
house is 30
+If confidence of meeting quota is zero
+(we assume all houses are successful sales), then finishing house is 5
+If confidence of meeting quota is 0, then finishing house is 5
+If confidence of meeting quota is 0.001, then finishing house is 5
+If confidence of meeting quota is 0.01, then finishing house is 5
+If confidence of meeting quota is 0.05, then finishing house is 6.2
+If confidence of meeting quota is 0.1, then finishing house is 7.06
+If confidence of meeting quota is 0.5, then finishing house is 11.3
+If confidence of meeting quota is 0.9, then finishing house is 17.8
+If confidence of meeting quota is 0.95, then finishing house is 20.1
+If confidence of meeting quota is 0.99, then finishing house is 24.8
+If confidence of meeting quota is 0.999, then finishing house is 31.1
+If confidence of meeting quota is 1, then finishing house is 1.#J
+House for 5th (last) sale. Probability (%)
+5 0.01024
+6 0.04096
+7 0.096256
+8 0.17367
+9 0.26657
+10 0.3669
+11 0.46723
+12 0.56182
+13 0.64696
+14 0.72074
+15 0.78272
+16 0.83343
+17 0.874
+18 0.90583
+19 0.93039
+20 0.94905
+21 0.96304
+22 0.97342
+23 0.98103
+24 0.98655
+25 0.99053
+26 0.99337
+27 0.99539
+28 0.99681
+29 0.9978
+30 0.99849
+
+*/
+
+
+
+
+
+
=======================================
--- /dev/null
+++ /trunk/libs/math/example/negative_binomial_example2.cpp Wed Jul 29
07:28:12 2009
@@ -0,0 +1,182 @@
+// negative_binomial_example2.cpp
+
+// Copyright Paul A. Bristow 2007.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Simple example demonstrating use of the Negative Binomial Distribution.
+
+#include <boost/math/distributions/negative_binomial.hpp>
+ using boost::math::negative_binomial_distribution;
+ using boost::math::negative_binomial; // typedef
+
+// In a sequence of trials or events
+// (Bernoulli, independent, yes or no, succeed or fail)
+// with success_fraction probability p,
+// negative_binomial is the probability that k or fewer failures
+// preceed the r th trial's success.
+
+#include <iostream>
+using std::cout;
+using std::endl;
+using std::setprecision;
+using std::showpoint;
+using std::setw;
+using std::left;
+using std::right;
+#include <limits>
+using std::numeric_limits;
+
+int main()
+{
+ cout << "Negative_binomial distribution - simple example 2" << endl;
+ // Construct a negative binomial distribution with:
+ // 8 successes (r), success fraction (p) 0.25 = 25% or 1 in 4 successes.
+ negative_binomial mynbdist(8, 0.25); // Shorter method using typedef.
+
+ // Display (to check) properties of the distribution just constructed.
+ cout << "mean(mynbdist) = " << mean(mynbdist) << endl; // 24
+ cout << "mynbdist.successes() = " << mynbdist.successes() << endl; // 8
+ // r th successful trial, after k failures, is r + k th trial.
+ cout << "mynbdist.success_fraction() = " << mynbdist.success_fraction()
<< endl;
+ // success_fraction = failures/successes or k/r = 0.25 or 25%.
+ cout << "mynbdist.percent success = " << mynbdist.success_fraction() *
100 << "%" << endl;
+ // Show as % too.
+ // Show some cumulative distribution function values for failures k = 2
and 8
+ cout << "cdf(mynbdist, 2.) = " << cdf(mynbdist, 2.) << endl; //
0.000415802001953125
+ cout << "cdf(mynbdist, 8.) = " << cdf(mynbdist, 8.) << endl; //
0.027129956288263202
+ cout << "cdf(complement(mynbdist, 8.)) = " << cdf(complement(mynbdist,
8.)) << endl; // 0.9728700437117368
+ // Check that cdf plus its complement is unity.
+ cout << "cdf + complement = " << cdf(mynbdist, 8.) +
cdf(complement(mynbdist, 8.)) << endl; // 1
+ // Note: No complement for pdf!
+
+ // Compare cdf with sum of pdfs.
+ double sum = 0.; // Calculate the sum of all the pdfs,
+ int k = 20; // for 20 failures
+ for(signed i = 0; i <= k; ++i)
+ {
+ sum += pdf(mynbdist, double(i));
+ }
+ // Compare with the cdf
+ double cdf8 = cdf(mynbdist, static_cast<double>(k));
+ double diff = sum - cdf8; // Expect the diference to be very small.
+ cout << setprecision(17) << "Sum pdfs = " << sum << ' ' // sum =
0.40025683281803698
+ << ", cdf = " << cdf(mynbdist, static_cast<double>(k)) // cdf =
0.40025683281803687
+ << ", difference = " // difference = 0.50000000000000000
+ << setprecision(1) << diff/ (std::numeric_limits<double>::epsilon() *
sum)
+ << " in epsilon units." << endl;
+
+ // Note: Use boost::math::tools::epsilon rather than std::numeric_limits
+ // to cover RealTypes that do not specialize numeric_limits.
+
+//[neg_binomial_example2
+
+ // Print a table of values that can be used to plot
+ // using Excel, or some other superior graphical display tool.
+
+ cout.precision(17); // Use max_digits10 precision, the maximum available
for a reference table.
+ cout << showpoint << endl; // include trailing zeros.
+ // This is a maximum possible precision for the type (here double) to
suit a reference table.
+ int maxk = static_cast<int>(2. * mynbdist.successes() /
mynbdist.success_fraction());
+ // This maxk shows most of the range of interest, probability about
0.0001 to 0.999.
+ cout << "\n"" k pdf cdf""\n" << endl;
+ for (int k = 0; k < maxk; k++)
+ {
+ cout << right << setprecision(17) << showpoint
+ << right << setw(3) << k << ", "
+ << left << setw(25) << pdf(mynbdist, static_cast<double>(k))
+ << left << setw(25) << cdf(mynbdist, static_cast<double>(k))
+ << endl;
+ }
+ cout << endl;
+//] [/ neg_binomial_example2]
+ return 0;
+} // int main()
+
+/*
+
+Output is:
+
+negative_binomial distribution - simple example 2
+mean(mynbdist) = 24
+mynbdist.successes() = 8
+mynbdist.success_fraction() = 0.25
+mynbdist.percent success = 25%
+cdf(mynbdist, 2.) = 0.000415802001953125
+cdf(mynbdist, 8.) = 0.027129956288263202
+cdf(complement(mynbdist, 8.)) = 0.9728700437117368
+cdf + complement = 1
+Sum pdfs = 0.40025683281803692 , cdf = 0.40025683281803687, difference =
0.25 in epsilon units.
+
+//[neg_binomial_example2_1
+ k pdf cdf
+ 0, 1.5258789062500000e-005 1.5258789062500003e-005
+ 1, 9.1552734375000000e-005 0.00010681152343750000
+ 2, 0.00030899047851562522 0.00041580200195312500
+ 3, 0.00077247619628906272 0.0011882781982421875
+ 4, 0.0015932321548461918 0.0027815103530883789
+ 5, 0.0028678178787231476 0.0056493282318115234
+ 6, 0.0046602040529251142 0.010309532284736633
+ 7, 0.0069903060793876605 0.017299838364124298
+ 8, 0.0098301179241389001 0.027129956288263202
+ 9, 0.013106823898851871 0.040236780187115073
+ 10, 0.016711200471036140 0.056947980658151209
+ 11, 0.020509200578089786 0.077457181236241013
+ 12, 0.024354675686481652 0.10181185692272265
+ 13, 0.028101548869017230 0.12991340579173993
+ 14, 0.031614242477644432 0.16152764826938440
+ 15, 0.034775666725408917 0.19630331499479325
+ 16, 0.037492515688331451 0.23379583068312471
+ 17, 0.039697957787645101 0.27349378847076977
+ 18, 0.041352039362130305 0.31484582783290005
+ 19, 0.042440250924291580 0.35728607875719176
+ 20, 0.042970754060845245 0.40025683281803687
+ 21, 0.042970754060845225 0.44322758687888220
+ 22, 0.042482450037426581 0.48571003691630876
+ 23, 0.041558918514873783 0.52726895543118257
+ 24, 0.040260202311284021 0.56752915774246648
+ 25, 0.038649794218832620 0.60617895196129912
+ 26, 0.036791631035234917 0.64297058299653398
+ 27, 0.034747651533277427 0.67771823452981139
+ 28, 0.032575923312447595 0.71029415784225891
+ 29, 0.030329307911589130 0.74062346575384819
+ 30, 0.028054609818219924 0.76867807557206813
+ 31, 0.025792141284492545 0.79447021685656061
+ 32, 0.023575629142856460 0.81804584599941710
+ 33, 0.021432390129869489 0.83947823612928651
+ 34, 0.019383705779220189 0.85886194190850684
+ 35, 0.017445335201298231 0.87630727710980494
+ 36, 0.015628112784496322 0.89193538989430121
+ 37, 0.013938587078064250 0.90587397697236549
+ 38, 0.012379666154859701 0.91825364312722524
+ 39, 0.010951243136991251 0.92920488626421649
+ 40, 0.0096507830144735539 0.93885566927869002
+ 41, 0.0084738582566109364 0.94732952753530097
+ 42, 0.0074146259745345548 0.95474415350983555
+ 43, 0.0064662435824429246 0.96121039709227851
+ 44, 0.0056212231142827853 0.96683162020656122
+ 45, 0.0048717266990450708 0.97170334690560634
+ 46, 0.0042098073105878630 0.97591315421619418
+ 47, 0.0036275999165703964 0.97954075413276465
+ 48, 0.0031174686783026818 0.98265822281106729
+ 49, 0.0026721160099737302 0.98533033882104104
+ 50, 0.0022846591885275322 0.98761499800956853
+ 51, 0.0019486798960970148 0.98956367790566557
+ 52, 0.0016582516423517923 0.99122192954801736
+ 53, 0.0014079495076571762 0.99262987905567457
+ 54, 0.0011928461106539983 0.99382272516632852
+ 55, 0.0010084971662802015 0.99483122233260868
+ 56, 0.00085091948404891532 0.99568214181665760
+ 57, 0.00071656377604119542 0.99639870559269883
+ 58, 0.00060228420831048650 0.99700098980100937
+ 59, 0.00050530624256557675 0.99750629604357488
+ 60, 0.00042319397814867202 0.99792949002172360
+ 61, 0.00035381791615708398 0.99828330793788067
+ 62, 0.00029532382517950324 0.99857863176306016
+ 63, 0.00024610318764958566 0.99882473495070978
+//] [neg_binomial_example2_1 end of Quickbook]
+
+*/
=======================================
--- /dev/null
+++ /trunk/libs/math/example/normal_misc_examples.cpp Wed Jul 29 07:28:12
2009
@@ -0,0 +1,509 @@
+// negative_binomial_example3.cpp
+
+// Copyright Paul A. Bristow 2007.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Example of using normal distribution.
+
+// Note that this file contains Quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+//[normal_basic1
+/*`
+First we need some includes to access the normal distribution
+(and some std output of course).
+*/
+
+#include <boost/math/distributions/normal.hpp> // for normal_distribution
+ using boost::math::normal; // typedef provides default type is double.
+
+#include <iostream>
+ using std::cout; using std::endl; using std::left; using std::showpoint;
using std::noshowpoint;
+#include <iomanip>
+ using std::setw; using std::setprecision;
+#include <limits>
+ using std::numeric_limits;
+
+int main()
+{
+ cout << "Example: Normal distribution, Miscellaneous Applications.";
+
+ try
+ {
+ { // Traditional tables and values.
+/*`Let's start by printing some traditional tables.
+*/
+ double step = 1.; // in z
+ double range = 4; // min and max z = -range to +range.
+ int precision = 17; // traditional tables are only computed to much
lower precision.
+
+ // Construct a standard normal distribution s
+ normal s; // (default mean = zero, and standard deviation = unity)
+ cout << "Standard normal distribution, mean = "<< s.mean()
+ << ", standard deviation = " << s.standard_deviation() << endl;
+
+/*` First the probability distribution function (pdf).
+*/
+ cout << "Probability distribution function values" << endl;
+ cout << " z " " pdf " << endl;
+ cout.precision(5);
+ for (double z = -range; z < range + step; z += step)
+ {
+ cout << left << setprecision(3) << setw(6) << z << " "
+ << setprecision(precision) << setw(12) << pdf(s, z) << endl;
+ }
+ cout.precision(6); // default
+ /*`And the area under the normal curve from -[infin] up to z,
+ the cumulative distribution function (cdf).
+*/
+ // For a standard normal distribution
+ cout << "Standard normal mean = "<< s.mean()
+ << ", standard deviation = " << s.standard_deviation() << endl;
+ cout << "Integral (area under the curve) from - infinity up to z "
<< endl;
+ cout << " z " " cdf " << endl;
+ for (double z = -range; z < range + step; z += step)
+ {
+ cout << left << setprecision(3) << setw(6) << z << " "
+ << setprecision(precision) << setw(12) << cdf(s, z) << endl;
+ }
+ cout.precision(6); // default
+
+/*`And all this you can do with a nanoscopic amount of work compared to
+the team of *human computers* toiling with Milton Abramovitz and Irene
Stegen
+at the US National Bureau of Standards (now [@http://www.nist.gov NIST]).
+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 [@http://dlmf.nist.gov Digital Library
of Mathematical Functions]).
+
+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
+*/
+ double z = 2.;
+ cout << "Area for z = " << z << " is " << cdf(s, z) << endl; // to get
the area for z.
+/*`
+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
+*/
+ cout << "95% of area has a z below " << quantile(s, 0.95) << endl;
+ // 95% of area has a z below 1.64485
+/*`and a two-sided test (a comparison between two levels, rather than a
one-sided test)
+
+*/
+ cout << "95% of area has a z between " << quantile(s, 0.975)
+ << " and " << -quantile(s, 0.975) << endl;
+ // 95% of area has a z between 1.95996 and -1.95996
+/*`
+
+First, define a table of significance levels: these are the probabilities
+that the true occurrence frequency lies outside the calculated interval.
+
+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,
+*/
+ double alpha1 = cdf(s, -1) * 2; // 0.3173105078629142
+ cout << setprecision(17) << "Significance level for z == 1 is " <<
alpha1 << endl;
+/*`
+ and place in our array of favorite alpha values.
+*/
+ double alpha[] = {0.3173105078629142, // z for 1 standard deviation.
+ 0.20, 0.1, 0.05, 0.01, 0.001, 0.0001, 0.00001 };
+/*`
+
+Confidence value as % is (1 - alpha) * 100 (so alpha 0.05 == 95%
confidence)
+that the true occurrence frequency lies *inside* the calculated interval.
+
+*/
+ cout << "level of significance (alpha)" << setprecision(4) << endl;
+ cout << "2-sided 1 -sided z(alpha) " << endl;
+ for (int i = 0; i < sizeof(alpha)/sizeof(alpha[0]); ++i)
+ {
+ cout << setw(15) << alpha[i] << setw(15) << alpha[i] /2 << setw(10)
<< quantile(complement(s, alpha[i]/2)) << endl;
+ // Use quantile(complement(s, alpha[i]/2)) to avoid potential loss
of accuracy from quantile(s, 1 - alpha[i]/2)
+ }
+ cout << endl;
+
+/*`Notice the distinction between one-sided (also called one-tailed)
+where we are using a > *or* < test (and not both)
+and considering the area of the tail (integral) from z up to +[infin],
+and a two-sided test where we are using two > *and* < tests, and thus
considering two tails,
+from -[infin] up to z low and z high up to +[infin].
+
+So the 2-sided values alpha[i] are calculated using alpha[i]/2.
+
+If we consider a simple example of alpha = 0.05, then for a two-sided test,
+the lower tail area from -[infin] 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,
+
+*/
+//] [/[normal_basic1]
+
+//[normal_basic2
+
+/*`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.
+
+*/
+ cout.precision(3);
+ cout << showpoint << "cdf(s, s.standard_deviation()) = "
+ << cdf(s, s.standard_deviation()) << endl; // from -infinity to 1 sd
+ cout << "cdf(complement(s, s.standard_deviation())) = "
+ << cdf(complement(s, s.standard_deviation())) << endl;
+ cout << "Fraction 1 standard deviation within either side of mean is "
+ << 1 - cdf(complement(s, s.standard_deviation())) * 2 << endl;
+ cout << "Fraction 2 standard deviations within either side of mean is "
+ << 1 - cdf(complement(s, 2 * s.standard_deviation())) * 2 << endl;
+ cout << "Fraction 3 standard deviations within either side of mean is "
+ << 1 - cdf(complement(s, 3 * s.standard_deviation())) * 2 << endl;
+
+/*`
+To a useful precision, the 1, 2 & 3 percentages are 68, 95 and 99.7,
+and these are worth memorising as useful 'rules of thumb', as, for
example, in
+[@http://en.wikipedia.org/wiki/Standard_deviation standard deviation]:
+
+[pre
+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
+]
+
+We could of course get some really accurate values for these
+[@http://en.wikipedia.org/wiki/Confidence_interval confidence intervals]
+by using cout.precision(15);
+
+[pre
+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
+]
+
+But before you get too excited about this impressive precision,
+don't forget that the *confidence intervals of the standard deviation* are
surprisingly wide,
+especially if you have estimated the standard deviation from only a few
measurements.
+*/
+//] [/[normal_basic2]
+
+
+//[normal_bulbs_example1
+/*`
+Examples from K. Krishnamoorthy, Handbook of Statistical Distributions
with Applications,
+ISBN 1 58488 635 8, page 125... implemented using the Math Toolkit library.
+
+A few very simple examples are shown here:
+*/
+// K. Krishnamoorthy, Handbook of Statistical Distributions with
Applications,
+ // ISBN 1 58488 635 8, page 125, example 10.3.5
+/*`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 /bulbs/ with these values:
+*/
+ double mean_life = 1100.;
+ double life_standard_deviation = 100.;
+ normal bulbs(mean_life, life_standard_deviation);
+ double expected_life = 1000.;
+
+/*`The we can use the Cumulative distribution function to predict fractions
+(or percentages, if * 100) that will last various lifetimes.
+*/
+ cout << "Fraction of bulbs that will last at best (<=) " // P(X <=
1000)
+ << expected_life << " is "<< cdf(bulbs, expected_life) << endl;
+ cout << "Fraction of bulbs that will last at least (>) " // P(X > 1000)
+ << expected_life << " is "<< cdf(complement(bulbs, expected_life))
<< endl;
+ double min_life = 900;
+ double max_life = 1200;
+ cout << "Fraction of bulbs that will last between "
+ << min_life << " and " << max_life << " is "
+ << cdf(bulbs, max_life) // P(X <= 1200)
+ - cdf(bulbs, min_life) << endl; // P(X <= 900)
+/*`
+[note 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.]
+*/
+//] [/normal_bulbs_example1 Quickbook end]
+ }
+ {
+ // K. Krishnamoorthy, Handbook of Statistical Distributions with
Applications,
+ // ISBN 1 58488 635 8, page 125, Example 10.3.6
+
+//[normal_bulbs_example3
+/*`Weekly demand for 5 lb sacks of onions at a store is normally
distributed with mean 140 sacks and standard deviation 10.
+*/
+ double mean = 140.; // sacks per week.
+ double standard_deviation = 10;
+ normal sacks(mean, standard_deviation);
+
+ double stock = 160.; // per week.
+ cout << "Percentage of weeks overstocked "
+ << cdf(sacks, stock) * 100. << endl; // P(X <=160)
+ // Percentage of weeks overstocked 97.7
+
+/*`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.
+*/
+ double stock_95 = quantile(sacks, 0.95);
+ cout << "Store should stock " << int(stock_95) << " sacks to meet 95% of
demands." << endl;
+/*`And it is easy to estimate how to meet 80% of demand, and waste even
less.
+*/
+ double stock_80 = quantile(sacks, 0.80);
+ cout << "Store should stock " << int(stock_80) << " sacks to meet 8 out
of 10 demands." << endl;
+//] [/normal_bulbs_example3 Quickbook end]
+ }
+ { // K. Krishnamoorthy, Handbook of Statistical Distributions with
Applications,
+ // ISBN 1 58488 635 8, page 125, Example 10.3.7
+
+//[normal_bulbs_example4
+
+/*`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.
+*/
+
+double mean = 3.; // kg
+double standard_deviation = 0.1; // kg
+normal packs(mean, standard_deviation);
+
+double max_weight = 3.1; // kg
+cout << "Percentage of packs > " << max_weight << " is "
+<< cdf(complement(packs, max_weight)) << endl; // P(X > 3.1)
+
+double under_weight = 2.9;
+cout <<"fraction of packs <= " << under_weight << " with a mean of " <<
mean
+ << " is " << cdf(complement(packs, under_weight)) << endl;
+// fraction of packs <= 2.9 with a mean of 3 is 0.841345
+// This is 0.84 - more than the target 0.95
+// Want 95% to be over this weight, so what should we set the mean weight
to be?
+// KK StatCalc says:
+double over_mean = 3.0664;
+normal xpacks(over_mean, standard_deviation);
+cout << "fraction of packs >= " << under_weight
+<< " with a mean of " << xpacks.mean()
+ << " is " << cdf(complement(xpacks, under_weight)) << endl;
+// fraction of packs >= 2.9 with a mean of 3.06449 is 0.950005
+double under_fraction = 0.05; // so 95% are above the minimum weight mean
- sd = 2.9
+double low_limit = standard_deviation;
+double offset = mean - low_limit - quantile(packs, under_fraction);
+double nominal_mean = mean + offset;
+
+normal nominal_packs(nominal_mean, standard_deviation);
+cout << "Setting the packer to " << nominal_mean << " will mean that "
+ << "fraction of packs >= " << under_weight
+ << " is " << cdf(complement(nominal_packs, under_weight)) << endl;
+
+/*`
+Setting the packer to 3.06449 will mean that fraction of packs >= 2.9 is
0.95.
+
+Setting the packer to 3.13263 will mean that fraction of packs >= 2.9 is
0.99,
+but will more than double the mean loss from 0.0644 to 0.133.
+
+Alternatively, we could invest in a better (more precise) packer with a
lower standard deviation.
+
+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
+*/
+double p = 0.05; // wanted p th quantile.
+cout << "Quantile of " << p << " = " << quantile(packs, p)
+ << ", mean = " << packs.mean() << ", sd = " <<
packs.standard_deviation() << endl; //
+/*`
+Quantile of 0.05 = 2.83551, mean = 3, sd = 0.1
+
+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.
+
+Let's start by guessing that it (now 0.1) needs to be halved, to a
standard deviation of 0.05
+*/
+normal pack05(mean, 0.05);
+cout << "Quantile of " << p << " = " << quantile(pack05, p)
+ << ", mean = " << pack05.mean() << ", sd = " <<
pack05.standard_deviation() << endl;
+
+cout <<"Fraction of packs >= " << under_weight << " with a mean of " <<
mean
+ << " and standard deviation of " << pack05.standard_deviation()
+ << " is " << cdf(complement(pack05, under_weight)) << endl;
+//
+/*`
+Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.05
is 0.9772
+
+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
+*/
+
+normal pack06(mean, 0.06);
+cout << "Quantile of " << p << " = " << quantile(pack06, p)
+ << ", mean = " << pack06.mean() << ", sd = " <<
pack06.standard_deviation() << endl;
+
+cout <<"Fraction of packs >= " << under_weight << " with a mean of " <<
mean
+ << " and standard deviation of " << pack06.standard_deviation()
+ << " is " << cdf(complement(pack06, under_weight)) << endl;
+/*`
+Fraction of packs >= 2.9 with a mean of 3 and standard deviation of 0.06
is 0.9522
+
+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
+[link math_toolkit.toolkit.internals1.roots2 Root Finding Without
Derivatives].
+
+But in this normal distribution case, we could be even smarter and make a
direct calculation.
+*/
+
+normal s; // For standard normal distribution,
+double sd = 0.1;
+double x = 2.9; // Our required limit.
+// then probability p = N((x - mean) / sd)
+// So if we want to find the standard deviation that would be required to
meet this limit,
+// so that the p th quantile is located at x,
+// in this case the 0.95 (95%) quantile at 2.9 kg pack weight, when the
mean is 3 kg.
+
+double prob = pdf(s, (x - mean) / sd);
+double qp = quantile(s, 0.95);
+cout << "prob = " << prob << ", quantile(p) " << qp << endl; // p =
0.241971, quantile(p) 1.64485
+// Rearranging, we can directly calculate the required standard deviation:
+double sd95 = abs((x - mean)) / qp;
+
+cout << "If we want the "<< p << " th quantile to be located at "
+ << x << ", would need a standard deviation of " << sd95 << endl;
+
+normal pack95(mean, sd95); // Distribution of the 'ideal better' packer.
+cout <<"Fraction of packs >= " << under_weight << " with a mean of " <<
mean
+ << " and standard deviation of " << pack95.standard_deviation()
+ << " is " << cdf(complement(pack95, under_weight)) << endl;
+
+// Fraction of packs >= 2.9 with a mean of 3 and standard deviation of
0.0608 is 0.95
+
+/*`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.
+*/
+
+//] [/normal_bulbs_example4 Quickbook end]
+ }
+
+ { // K. Krishnamoorthy, Handbook of Statistical Distributions with
Applications,
+ // ISBN 1 58488 635 8, page 125, example 10.3.8
+//[normal_bulbs_example5
+/*`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.
+*/
+
+ normal bolts(3.95, 0.1);
+ double top = 4.1;
+ double bottom = 3.9;
+
+cout << "Fraction long enough [ P(X <= " << top << ") ] is " << cdf(bolts,
top) << endl;
+cout << "Fraction too short [ P(X <= " << bottom << ") ] is " <<
cdf(bolts, bottom) << endl;
+cout << "Fraction OK -between " << bottom << " and " << top
+ << "[ P(X <= " << top << ") - P(X<= " << bottom << " ) ] is "
+ << cdf(bolts, top) - cdf(bolts, bottom) << endl;
+
+cout << "Fraction too long [ P(X > " << top << ") ] is "
+ << cdf(complement(bolts, top)) << endl;
+
+cout << "95% of bolts are shorter than " << quantile(bolts, 0.95) << endl;
+
+//] [/normal_bulbs_example5 Quickbook end]
+ }
+ }
+ catch(const std::exception& e)
+ { // Always useful to include try & catch blocks because default policies
+ // are to throw exceptions on arguments that cause errors like
underflow, overflow.
+ // Lacking try & catch blocks, the program will abort without a
message below,
+ // which may give some helpful clues as to the cause of the exception.
+ std::cout <<
+ "\n""Message from thrown exception was:\n " << e.what() <<
std::endl;
+ }
+ return 0;
+} // int main()
+
+
+/*
+
+Output is:
+
+Autorun
"i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\normal_misc_examples.exe"
+Example: Normal distribution, Miscellaneous Applications.Standard normal
distribution, mean = 0, standard deviation = 1
+Probability distribution function values
+ z pdf
+-4 0.00013383022576488537
+-3 0.0044318484119380075
+-2 0.053990966513188063
+-1 0.24197072451914337
+0 0.3989422804014327
+1 0.24197072451914337
+2 0.053990966513188063
+3 0.0044318484119380075
+4 0.00013383022576488537
+Standard normal mean = 0, standard deviation = 1
+Integral (area under the curve) from - infinity up to z
+ z cdf
+-4 3.1671241833119979e-005
+-3 0.0013498980316300959
+-2 0.022750131948179219
+-1 0.1586552539314571
+0 0.5
+1 0.84134474606854293
+2 0.97724986805182079
+3 0.9986501019683699
+4 0.99996832875816688
+Area for z = 2 is 0.97725
+95% of area has a z below 1.64485
+95% of area has a z between 1.95996 and -1.95996
+Significance level for z == 1 is 0.3173105078629142
+level of significance (alpha)
+2-sided 1 -sided z(alpha)
+0.3173 0.1587 1
+0.2 0.1 1.282
+0.1 0.05 1.645
+0.05 0.025 1.96
+0.01 0.005 2.576
+0.001 0.0005 3.291
+0.0001 5e-005 3.891
+1e-005 5e-006 4.417
+cdf(s, s.standard_deviation()) = 0.841
+cdf(complement(s, s.standard_deviation())) = 0.159
+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
+Fraction of bulbs that will last at best (<=) 1.00e+003 is 0.159
+Fraction of bulbs that will last at least (>) 1.00e+003 is 0.841
+Fraction of bulbs that will last between 900. and 1.20e+003 is 0.819
+Percentage of weeks overstocked 97.7
+Store should stock 156 sacks to meet 95% of demands.
+Store should stock 148 sacks to meet 8 out of 10 demands.
+Percentage of packs > 3.10 is 0.159
+fraction of packs <= 2.90 with a mean of 3.00 is 0.841
+fraction of packs >= 2.90 with a mean of 3.07 is 0.952
+Setting the packer to 3.06 will mean that fraction of packs >= 2.90 is
0.950
+Quantile of 0.0500 = 2.84, mean = 3.00, sd = 0.100
+Quantile of 0.0500 = 2.92, mean = 3.00, sd = 0.0500
+Fraction of packs >= 2.90 with a mean of 3.00 and standard deviation of
0.0500 is 0.977
+Quantile of 0.0500 = 2.90, mean = 3.00, sd = 0.0600
+Fraction of packs >= 2.90 with a mean of 3.00 and standard deviation of
0.0600 is 0.952
+prob = 0.242, quantile(p) 1.64
+If we want the 0.0500 th quantile to be located at 2.90, would need a
standard deviation of 0.0608
+Fraction of packs >= 2.90 with a mean of 3.00 and standard deviation of
0.0608 is 0.950
+Fraction long enough [ P(X <= 4.10) ] is 0.933
+Fraction too short [ P(X <= 3.90) ] is 0.309
+Fraction OK -between 3.90 and 4.10[ P(X <= 4.10) - P(X<= 3.90 ) ] is 0.625
+Fraction too long [ P(X > 4.10) ] is 0.0668
+95% of bolts are shorter than 4.11
+
+*/
+
+
+
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_eg_1.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,35 @@
+// Copyright John Maddock 2007.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include <iostream>
+
+//[policy_eg_1
+
+#include <boost/math/special_functions/gamma.hpp>
+//
+// Define the policy to use:
+//
+using namespace boost::math::policies;
+typedef policy<
+ domain_error<errno_on_error>,
+ pole_error<errno_on_error>,
+ overflow_error<errno_on_error>,
+ evaluation_error<errno_on_error>
+> c_policy;
+//
+// Now use the policy when calling tgamma:
+//
+int main()
+{
+ errno = 0;
+ std::cout << "Result of tgamma(30000) is: "
+ << boost::math::tgamma(30000, c_policy()) << std::endl;
+ std::cout << "errno = " << errno << std::endl;
+ std::cout << "Result of tgamma(-10) is: "
+ << boost::math::tgamma(-10, c_policy()) << std::endl;
+ std::cout << "errno = " << errno << std::endl;
+}
+
+//]
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_eg_10.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,171 @@
+// Copyright John Maddock 2007.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Note that this file contains quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+//[policy_eg_10
+
+/*`
+
+To understand how the rounding policies for
+the discrete distributions can be used, we'll
+use the 50-sample binomial distribution with a
+success fraction of 0.5 once again, and calculate
+all the possible quantiles at 0.05 and 0.95.
+
+Begin by including the needed headers:
+
+*/
+
+#include <iostream>
+#include <boost/math/distributions/binomial.hpp>
+
+/*`
+
+Next we'll bring the needed declarations into scope, and
+define distribution types for all the available rounding policies:
+
+*/
+
+using namespace boost::math::policies;
+using namespace boost::math;
+
+typedef binomial_distribution<
+ double,
+ policy<discrete_quantile<integer_round_outwards> > >
+ binom_round_outwards;
+
+typedef binomial_distribution<
+ double,
+ policy<discrete_quantile<integer_round_inwards> > >
+ binom_round_inwards;
+
+typedef binomial_distribution<
+ double,
+ policy<discrete_quantile<integer_round_down> > >
+ binom_round_down;
+
+typedef binomial_distribution<
+ double,
+ policy<discrete_quantile<integer_round_up> > >
+ binom_round_up;
+
+typedef binomial_distribution<
+ double,
+ policy<discrete_quantile<integer_round_nearest> > >
+ binom_round_nearest;
+
+typedef binomial_distribution<
+ double,
+ policy<discrete_quantile<real> > >
+ binom_real_quantile;
+
+/*`
+
+Now let's set to work calling those quantiles:
+
+*/
+
+int main()
+{
+ std::cout <<
+ "Testing rounding policies for a 50 sample binomial distribution,\n"
+ "with a success fraction of 0.5.\n\n"
+ "Lower quantiles are calculated at p = 0.05\n\n"
+ "Upper quantiles at p = 0.95.\n\n";
+
+ std::cout << std::setw(25) << std::right
+ << "Policy"<< std::setw(18) << std::right
+ << "Lower Quantile" << std::setw(18) << std::right
+ << "Upper Quantile" << std::endl;
+
+ // Test integer_round_outwards:
+ std::cout << std::setw(25) << std::right
+ << "integer_round_outwards"
+ << std::setw(18) << std::right
+ << quantile(binom_round_outwards(50, 0.5), 0.05)
+ << std::setw(18) << std::right
+ << quantile(binom_round_outwards(50, 0.5), 0.95)
+ << std::endl;
+
+ // Test integer_round_inwards:
+ std::cout << std::setw(25) << std::right
+ << "integer_round_inwards"
+ << std::setw(18) << std::right
+ << quantile(binom_round_inwards(50, 0.5), 0.05)
+ << std::setw(18) << std::right
+ << quantile(binom_round_inwards(50, 0.5), 0.95)
+ << std::endl;
+
+ // Test integer_round_down:
+ std::cout << std::setw(25) << std::right
+ << "integer_round_down"
+ << std::setw(18) << std::right
+ << quantile(binom_round_down(50, 0.5), 0.05)
+ << std::setw(18) << std::right
+ << quantile(binom_round_down(50, 0.5), 0.95)
+ << std::endl;
+
+ // Test integer_round_up:
+ std::cout << std::setw(25) << std::right
+ << "integer_round_up"
+ << std::setw(18) << std::right
+ << quantile(binom_round_up(50, 0.5), 0.05)
+ << std::setw(18) << std::right
+ << quantile(binom_round_up(50, 0.5), 0.95)
+ << std::endl;
+
+ // Test integer_round_nearest:
+ std::cout << std::setw(25) << std::right
+ << "integer_round_nearest"
+ << std::setw(18) << std::right
+ << quantile(binom_round_nearest(50, 0.5), 0.05)
+ << std::setw(18) << std::right
+ << quantile(binom_round_nearest(50, 0.5), 0.95)
+ << std::endl;
+
+ // Test real:
+ std::cout << std::setw(25) << std::right
+ << "real"
+ << std::setw(18) << std::right
+ << quantile(binom_real_quantile(50, 0.5), 0.05)
+ << std::setw(18) << std::right
+ << quantile(binom_real_quantile(50, 0.5), 0.95)
+ << std::endl;
+}
+
+/*`
+
+Which produces the program output:
+
+[pre
+Testing rounding policies for a 50 sample binomial distribution,
+with a success fraction of 0.5.
+
+Lower quantiles are calculated at p = 0.05
+
+Upper quantiles at p = 0.95.
+
+Testing rounding policies for a 50 sample binomial distribution,
+with a success fraction of 0.5.
+
+Lower quantiles are calculated at p = 0.05
+
+Upper quantiles at p = 0.95.
+
+ Policy Lower Quantile Upper Quantile
+ integer_round_outwards 18 31
+ integer_round_inwards 19 30
+ integer_round_down 18 30
+ integer_round_up 19 31
+ integer_round_nearest 19 30
+ real 18.701 30.299
+]
+
+*/
+
+//] ends quickbook import
+
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_eg_2.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,44 @@
+// Copyright John Maddock 2007.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#include <iostream>
+
+//[policy_eg_2
+
+#include <boost/math/special_functions/gamma.hpp>
+
+int main()
+{
+ using namespace boost::math::policies;
+ errno = 0;
+ std::cout << "Result of tgamma(30000) is: "
+ << boost::math::tgamma(
+ 30000,
+ make_policy(
+ domain_error<errno_on_error>(),
+ pole_error<errno_on_error>(),
+ overflow_error<errno_on_error>(),
+ evaluation_error<errno_on_error>()
+ )
+ ) << std::endl;
+ // Check errno was set:
+ std::cout << "errno = " << errno << std::endl;
+ // and again with evaluation at a pole:
+ std::cout << "Result of tgamma(-10) is: "
+ << boost::math::tgamma(
+ -10,
+ make_policy(
+ domain_error<errno_on_error>(),
+ pole_error<errno_on_error>(),
+ overflow_error<errno_on_error>(),
+ evaluation_error<errno_on_error>()
+ )
+ ) << std::endl;
+ // Check errno was set:
+ std::cout << "errno = " << errno << std::endl;
+}
+
+//]
+
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_eg_3.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,42 @@
+// Copyright John Maddock 2007.
+// Copyright Paul A. Bristow 2007.
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+#ifdef _MSC_VER
+# pragma warning (disable : 4305) // 'initializing' : truncation
from 'long double' to 'const eval_type'
+# pragma warning (disable : 4244) // conversion from 'long double'
to 'const eval_type'
+#endif
+
+#include <iostream>
+
+//[policy_eg_3
+
+#include <boost/math/distributions/binomial.hpp>
+
+//
+// Begin by defining a policy type, that gives the
+// behaviour we want:
+//
+using namespace boost::math::policies;
+typedef policy<
+ promote_float<false>,
+ discrete_quantile<integer_round_nearest>
+> mypolicy;
+//
+// Then define a distribution that uses it:
+//
+typedef boost::math::binomial_distribution<float, mypolicy> mybinom;
+//
+// And now use it to get the quantile:
+//
+int main()
+{
+ std::cout << "quantile is: " <<
+ quantile(mybinom(200, 0.25), 0.05) << std::endl;
+}
+
+//]
+
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_eg_4.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,102 @@
+// Copyright John Maddock 2007.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Note that this file contains quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+#include <iostream>
+
+//[policy_eg_4
+
+/*`
+Suppose we want `C::foo()` to behave in a C-compatible way and set
+`::errno` on error rather than throwing any exceptions.
+
+We'll begin by including the needed header:
+*/
+
+#include <boost/math/special_functions.hpp>
+
+/*`
+
+Open up the "C" namespace that we'll use for our functions, and
+define the policy type we want: in this case one that sets
+::errno rather than throwing exceptions. Any policies we don't
+specify here will inherit the defaults:
+
+*/
+
+namespace C{
+
+using namespace boost::math::policies;
+
+typedef policy<
+ domain_error<errno_on_error>,
+ pole_error<errno_on_error>,
+ overflow_error<errno_on_error>,
+ evaluation_error<errno_on_error>
+> c_policy;
+
+/*`
+
+All we need do now is invoke the BOOST_MATH_DECLARE_SPECIAL_FUNCTIONS
+macro passing our policy type as the single argument:
+
+*/
+
+BOOST_MATH_DECLARE_SPECIAL_FUNCTIONS(c_policy)
+
+} // close namespace C
+
+/*`
+
+We now have a set of forwarding functions defined in namespace C
+that all look something like this:
+
+``
+template <class RealType>
+inline typename boost::math::tools::promote_args<RT>::type
+ tgamma(RT z)
+{
+ return boost::math::tgamma(z, c_policy());
+}
+``
+
+So that when we call `C::tgamma(z)` we really end up calling
+`boost::math::tgamma(z, C::c_policy())`:
+
+*/
+
+
+int main()
+{
+ errno = 0;
+ std::cout << "Result of tgamma(30000) is: "
+ << C::tgamma(30000) << std::endl;
+ std::cout << "errno = " << errno << std::endl;
+ std::cout << "Result of tgamma(-10) is: "
+ << C::tgamma(-10) << std::endl;
+ std::cout << "errno = " << errno << std::endl;
+}
+
+/*`
+
+Which outputs:
+
+[pre
+Result of C::tgamma(30000) is: 1.#INF
+errno = 34
+Result of C::tgamma(-10) is: 1.#QNAN
+errno = 33
+]
+
+This mechanism is particularly useful when we want to define a
+project-wide policy, and don't want to modify the Boost source
+or set - possibly fragile and easy to forget - project wide
+build macros.
+
+*/
+//]
+
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_eg_5.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,41 @@
+// Copyright John Maddock 2007.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Note that this file contains quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+#include <iostream>
+
+//[policy_eg_5
+
+#include <boost/math/special_functions.hpp>
+
+namespace {
+
+using namespace boost::math::policies;
+
+typedef policy<
+ domain_error<errno_on_error>,
+ pole_error<errno_on_error>,
+ overflow_error<errno_on_error>,
+ evaluation_error<errno_on_error>
+> c_policy;
+
+BOOST_MATH_DECLARE_SPECIAL_FUNCTIONS(c_policy)
+
+} // close unnamed namespace
+
+int main()
+{
+ errno = 0;
+ std::cout << "Result of tgamma(30000) is: "
+ << tgamma(30000) << std::endl;
+ std::cout << "errno = " << errno << std::endl;
+ std::cout << "Result of tgamma(-10) is: "
+ << tgamma(-10) << std::endl;
+ std::cout << "errno = " << errno << std::endl;
+}
+
+//]
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_eg_6.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,119 @@
+// Copyright John Maddock 2007.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Note that this file contains quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+#include <iostream>
+
+//[policy_eg_6
+
+/*`
+Suppose we want a set of distributions to behave as follows:
+
+* Return infinity on overflow, rather than throwing an exception.
+* Don't perform any promotion from double to long double internally.
+* Return the closest integer result from the quantiles of discrete
+distributions.
+
+We'll begin by including the needed header:
+*/
+
+#include <boost/math/distributions.hpp>
+
+/*`
+
+Open up an appropriate namespace for our distributions, and
+define the policy type we want. Any policies we don't
+specify here will inherit the defaults:
+
+*/
+
+namespace my_distributions{
+
+using namespace boost::math::policies;
+
+typedef policy<
+ // return infinity and set errno rather than throw:
+ overflow_error<errno_on_error>,
+ // Don't promote double -> long double internally:
+ promote_double<false>,
+ // Return the closest integer result for discrete quantiles:
+ discrete_quantile<integer_round_nearest>
+> my_policy;
+
+/*`
+
+All we need do now is invoke the BOOST_MATH_DECLARE_DISTRIBUTIONS
+macro passing the floating point and policy types as arguments:
+
+*/
+
+BOOST_MATH_DECLARE_DISTRIBUTIONS(double, my_policy)
+
+} // close namespace my_namespace
+
+/*`
+
+We now have a set of typedefs defined in namespace my_namespace
+that all look something like this:
+
+``
+typedef boost::math::normal_distribution<double, my_policy> normal;
+typedef boost::math::cauchy_distribution<double, my_policy> cauchy;
+typedef boost::math::gamma_distribution<double, my_policy> gamma;
+// etc
+``
+
+So that when we use my_namespace::normal we really end up using
+`boost::math::normal_distribution<double, my_policy>`:
+
+*/
+
+
+int main()
+{
+ //
+ // Start with something we know will overflow:
+ //
+ my_distributions::normal norm(10, 2);
+ errno = 0;
+ std::cout << "Result of quantile(norm, 0) is: "
+ << quantile(norm, 0) << std::endl;
+ std::cout << "errno = " << errno << std::endl;
+ errno = 0;
+ std::cout << "Result of quantile(norm, 1) is: "
+ << quantile(norm, 1) << std::endl;
+ std::cout << "errno = " << errno << std::endl;
+ //
+ // Now try a discrete distribution:
+ //
+ my_distributions::binomial binom(20, 0.25);
+ std::cout << "Result of quantile(binom, 0.05) is: "
+ << quantile(binom, 0.05) << std::endl;
+ std::cout << "Result of quantile(complement(binom, 0.05)) is: "
+ << quantile(complement(binom, 0.05)) << std::endl;
+}
+
+/*`
+
+Which outputs:
+
+[pre
+Result of quantile(norm, 0) is: -1.#INF
+errno = 34
+Result of quantile(norm, 1) is: 1.#INF
+errno = 34
+Result of quantile(binom, 0.05) is: 1
+Result of quantile(complement(binom, 0.05)) is: 8
+]
+
+This mechanism is particularly useful when we want to define a
+project-wide policy, and don't want to modify the Boost source
+or set - possibly fragile and easy to forget - project wide
+build macros.
+
+*/ //] ends quickbook imported section
+
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_eg_7.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,57 @@
+// Copyright John Maddock 2007.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Note that this file contains quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+#include <iostream>
+
+//[policy_eg_7
+
+#include <boost/math/distributions.hpp>
+
+namespace {
+
+using namespace boost::math::policies;
+
+typedef policy<
+ // return infinity and set errno rather than throw:
+ overflow_error<errno_on_error>,
+ // Don't promote double -> long double internally:
+ promote_double<false>,
+ // Return the closest integer result for discrete quantiles:
+ discrete_quantile<integer_round_nearest>
+> my_policy;
+
+BOOST_MATH_DECLARE_DISTRIBUTIONS(double, my_policy)
+
+} // close namespace my_namespace
+
+int main()
+{
+ //
+ // Start with something we know will overflow:
+ //
+ normal norm(10, 2);
+ errno = 0;
+ std::cout << "Result of quantile(norm, 0) is: "
+ << quantile(norm, 0) << std::endl;
+ std::cout << "errno = " << errno << std::endl;
+ errno = 0;
+ std::cout << "Result of quantile(norm, 1) is: "
+ << quantile(norm, 1) << std::endl;
+ std::cout << "errno = " << errno << std::endl;
+ //
+ // Now try a discrete distribution:
+ //
+ binomial binom(20, 0.25);
+ std::cout << "Result of quantile(binom, 0.05) is: "
+ << quantile(binom, 0.05) << std::endl;
+ std::cout << "Result of quantile(complement(binom, 0.05)) is: "
+ << quantile(complement(binom, 0.05)) << std::endl;
+}
+
+//] ends quickbook imported section
+
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_eg_8.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,137 @@
+// Copyright John Maddock 2007.
+// Copyright Paul a. Bristow 2007
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Note that this file contains quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+#ifdef _MSC_VER
+# pragma warning (disable : 4100) // 'unreferenced formal parameter
+#endif
+
+
+#include <iostream>
+
+//[policy_eg_8
+
+/*`
+
+Suppose we want our own user-defined error handlers rather than the
+any of the default ones supplied by the library to be used. If
+we set the policy for a specific type of error to `user_error`
+then the library will call a user-supplied error handler.
+These are forward declared, but not defined in
+boost/math/policies/error_handling.hpp like this:
+
+ namespace boost{ namespace math{ namespace policies{
+
+ template <class T>
+ T user_domain_error(const char* function, const char* message, const T&
val);
+ template <class T>
+ T user_pole_error(const char* function, const char* message, const T&
val);
+ template <class T>
+ T user_overflow_error(const char* function, const char* message, const
T& val);
+ template <class T>
+ T user_underflow_error(const char* function, const char* message, const
T& val);
+ template <class T>
+ T user_denorm_error(const char* function, const char* message, const T&
val);
+ template <class T>
+ T user_evaluation_error(const char* function, const char* message,
const T& val);
+ template <class T>
+ T user_indeterminate_result_error(const char* function, const char*
message, const T& val);
+
+ }}} // namespaces
+
+So out first job is to include the header we want to use, and then
+provide definitions for the user-defined error handlers we want to use:
+
+*/
+
+#include <iostream>
+#include <boost/math/special_functions.hpp>
+
+namespace boost{ namespace math{ namespace policies{
+
+template <class T>
+T user_domain_error(const char* function, const char* message, const T&
val)
+{
+ std::cerr << "Domain Error." << std::endl;
+ return std::numeric_limits<T>::quiet_NaN();
+}
+
+template <class T>
+T user_pole_error(const char* function, const char* message, const T& val)
+{
+ std::cerr << "Pole Error." << std::endl;
+ return std::numeric_limits<T>::quiet_NaN();
+}
+
+
+}}} // namespaces
+
+
+/*`
+
+Now we'll need to define a suitable policy that will call these handlers,
+and define some forwarding functions that make use of the policy:
+
+*/
+
+namespace{
+
+using namespace boost::math::policies;
+
+typedef policy<
+ domain_error<user_error>,
+ pole_error<user_error>
+> user_error_policy;
+
+BOOST_MATH_DECLARE_SPECIAL_FUNCTIONS(user_error_policy)
+
+} // close unnamed namespace
+
+/*`
+
+We now have a set of forwarding functions defined in an unnamed namespace
+that all look something like this:
+
+``
+template <class RealType>
+inline typename boost::math::tools::promote_args<RT>::type
+ tgamma(RT z)
+{
+ return boost::math::tgamma(z, user_error_policy());
+}
+``
+
+So that when we call `tgamma(z)` we really end up calling
+`boost::math::tgamma(z, user_error_policy())`, and any
+errors will get directed to our own error handlers:
+
+*/
+
+
+int main()
+{
+ std::cout << "Result of erf_inv(-10) is: "
+ << erf_inv(-10) << std::endl;
+ std::cout << "Result of tgamma(-10) is: "
+ << tgamma(-10) << std::endl;
+}
+
+/*`
+
+Which outputs:
+
+[pre
+Domain Error.
+Result of erf_inv(-10) is: 1.#QNAN
+Pole Error.
+Result of tgamma(-10) is: 1.#QNAN
+]
+*/
+
+//]
+
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_eg_9.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,316 @@
+// Copyright John Maddock 2007.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Note that this file contains quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+#include <iostream>
+
+//[policy_eg_9
+
+/*`
+
+The previous example was all well and good, but the custom error handlers
+didn't really do much of any use. In this example we'll implement all
+the custom handlers and show how the information provided to them can be
+used to generate nice formatted error messages.
+
+Each error handler has the general form:
+
+ template <class T>
+ T user_``['error_type]``(
+ const char* function,
+ const char* message,
+ const T& val);
+
+and accepts three arguments:
+
+[variablelist
+[[const char* function]
+ [The name of the function that raised the error, this string
+ contains one or more %1% format specifiers that should be
+ replaced by the name of type T.]]
+[[const char* message]
+ [A message associated with the error, normally this
+ contains a %1% format specifier that should be replaced with
+ the value of ['value]: however note that overflow and underflow messages
+ do not contain this %1% specifier (since the value of ['value] is
+ immaterial in these cases).]]
+[[const T& value]
+ [The value that caused the error: either an argument to the function
+ if this is a domain or pole error, the tentative result
+ if this is a denorm or evaluation error, or zero or infinity for
+ underflow or overflow errors.]]
+]
+
+As before we'll include the headers we need first:
+
+*/
+
+#include <iostream>
+#include <boost/math/special_functions.hpp>
+
+/*`
+
+Next we'll implement the error handlers for each type of error,
+starting with domain errors:
+
+*/
+
+namespace boost{ namespace math{ namespace policies{
+
+template <class T>
+T user_domain_error(const char* function, const char* message, const T&
val)
+{
+ /*`
+ We'll begin with a bit of defensive programming:
+ */
+ if(function == 0)
+ function = "Unknown function with arguments of type %1%";
+ if(message == 0)
+ message = "Cause unknown with bad argument %1%";
+ /*`
+ Next we'll format the name of the function with the name of type T:
+ */
+ std::string msg("Error in function ");
+ msg += (boost::format(function) % typeid(T).name()).str();
+ /*`
+ Then likewise format the error message with the value of parameter
/val/,
+ making sure we output all the digits of /val/:
+ */
+ msg += ": \n";
+ int prec = 2 + (std::numeric_limits<T>::digits * 30103UL) / 100000UL;
+ msg += (boost::format(message) %
boost::io::group(std::setprecision(prec), val)).str();
+ /*`
+ Now we just have to do something with the message, we could throw an
+ exception, but for the purposes of this example we'll just dump the
message
+ to std::cerr:
+ */
+ std::cerr << msg << std::endl;
+ /*`
+ Finally the only sensible value we can return from a domain error is a
NaN:
+ */
+ return std::numeric_limits<T>::quiet_NaN();
+}
+
+/*`
+Pole errors are essentially a special case of domain errors,
+so in this example we'll just return the result of a domain error:
+*/
+
+template <class T>
+T user_pole_error(const char* function, const char* message, const T& val)
+{
+ return user_domain_error(function, message, val);
+}
+
+/*`
+Overflow errors are very similar to domain errors, except that there's
+no %1% format specifier in the /message/ parameter:
+*/
+template <class T>
+T user_overflow_error(const char* function, const char* message, const T&
val)
+{
+ if(function == 0)
+ function = "Unknown function with arguments of type %1%";
+ if(message == 0)
+ message = "Result of function is too large to represent";
+
+ std::string msg("Error in function ");
+ msg += (boost::format(function) % typeid(T).name()).str();
+
+ msg += ": \n";
+ msg += message;
+
+ std::cerr << msg << std::endl;
+
+ // Value passed to the function is an infinity, just return it:
+ return val;
+}
+
+/*`
+Underflow errors are much the same as overflow:
+*/
+
+template <class T>
+T user_underflow_error(const char* function, const char* message, const T&
val)
+{
+ if(function == 0)
+ function = "Unknown function with arguments of type %1%";
+ if(message == 0)
+ message = "Result of function is too small to represent";
+
+ std::string msg("Error in function ");
+ msg += (boost::format(function) % typeid(T).name()).str();
+
+ msg += ": \n";
+ msg += message;
+
+ std::cerr << msg << std::endl;
+
+ // Value passed to the function is zero, just return it:
+ return val;
+}
+
+/*`
+Denormalised results are much the same as underflow:
+*/
+
+template <class T>
+T user_denorm_error(const char* function, const char* message, const T&
val)
+{
+ if(function == 0)
+ function = "Unknown function with arguments of type %1%";
+ if(message == 0)
+ message = "Result of function is denormalised";
+
+ std::string msg("Error in function ");
+ msg += (boost::format(function) % typeid(T).name()).str();
+
+ msg += ": \n";
+ msg += message;
+
+ std::cerr << msg << std::endl;
+
+ // Value passed to the function is denormalised, just return it:
+ return val;
+}
+
+/*`
+Which leaves us with evaluation errors, these occur when an internal
+error occurs that prevents the function being fully evaluated.
+The parameter /val/ contains the closest approximation to the result
+found so far:
+*/
+
+template <class T>
+T user_evaluation_error(const char* function, const char* message, const
T& val)
+{
+ if(function == 0)
+ function = "Unknown function with arguments of type %1%";
+ if(message == 0)
+ message = "An internal evaluation error occured with "
+ "the best value calculated so far of %1%";
+
+ std::string msg("Error in function ");
+ msg += (boost::format(function) % typeid(T).name()).str();
+
+ msg += ": \n";
+ int prec = 2 + (std::numeric_limits<T>::digits * 30103UL) / 100000UL;
+ msg += (boost::format(message) %
boost::io::group(std::setprecision(prec), val)).str();
+
+ std::cerr << msg << std::endl;
+
+ // What do we return here? This is generally a fatal error,
+ // that should never occur, just return a NaN for the purposes
+ // of the example:
+ return std::numeric_limits<T>::quiet_NaN();
+}
+
+}}} // namespaces
+
+
+/*`
+
+Now we'll need to define a suitable policy that will call these handlers,
+and define some forwarding functions that make use of the policy:
+
+*/
+
+namespace{
+
+using namespace boost::math::policies;
+
+typedef policy<
+ domain_error<user_error>,
+ pole_error<user_error>,
+ overflow_error<user_error>,
+ underflow_error<user_error>,
+ denorm_error<user_error>,
+ evaluation_error<user_error>
+> user_error_policy;
+
+BOOST_MATH_DECLARE_SPECIAL_FUNCTIONS(user_error_policy)
+
+} // close unnamed namespace
+
+/*`
+
+We now have a set of forwarding functions defined in an unnamed namespace
+that all look something like this:
+
+``
+template <class RealType>
+inline typename boost::math::tools::promote_args<RT>::type
+ tgamma(RT z)
+{
+ return boost::math::tgamma(z, user_error_policy());
+}
+``
+
+So that when we call `tgamma(z)` we really end up calling
+`boost::math::tgamma(z, user_error_policy())`, and any
+errors will get directed to our own error handlers:
+
+*/
+
+
+int main()
+{
+ // Raise a domain error:
+ std::cout << "Result of erf_inv(-10) is: "
+ << erf_inv(-10) << std::endl << std::endl;
+ // Raise a pole error:
+ std::cout << "Result of tgamma(-10) is: "
+ << tgamma(-10) << std::endl << std::endl;
+ // Raise an overflow error:
+ std::cout << "Result of tgamma(3000) is: "
+ << tgamma(3000) << std::endl << std::endl;
+ // Raise an underflow error:
+ std::cout << "Result of tgamma(-190.5) is: "
+ << tgamma(-190.5) << std::endl << std::endl;
+ // Unfortunately we can't predicably raise a denormalised
+ // result, nor can we raise an evaluation error in this example
+ // since these should never really occur!
+}
+
+/*`
+
+Which outputs:
+
+[pre
+Error in function boost::math::erf_inv<double>(double, double):
+Argument outside range \[-1, 1\] in inverse erf function (got p=-10).
+Result of erf_inv(-10) is: 1.#QNAN
+
+Error in function boost::math::tgamma<long double>(long double):
+Evaluation of tgamma at a negative integer -10.
+Result of tgamma(-10) is: 1.#QNAN
+
+Error in function boost::math::tgamma<long double>(long double):
+Result of tgamma is too large to represent.
+Error in function boost::math::tgamma<double>(double):
+Result of function is too large to represent
+Result of tgamma(3000) is: 1.#INF
+
+Error in function boost::math::tgamma<long double>(long double):
+Result of tgamma is too large to represent.
+Error in function boost::math::tgamma<long double>(long double):
+Result of tgamma is too small to represent.
+Result of tgamma(-190.5) is: 0
+]
+
+Notice how some of the calls result in an error handler being called more
+than once, or for more than one handler to be called: this is an artefact
+of the fact that many functions are implemented in terms of one or more
+sub-routines each of which may have it's own error handling. For example
+`tgamma(-190.5)` is implemented in terms of `tgamma(190.5)` - which
overflows -
+the reflection formula for `tgamma` then notices that it's dividing by
+infinity and underflows.
+
+*/
+
+//]
+
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_ref_snip1.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,44 @@
+// Copyright John Maddock 2007.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Note that this file contains quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+double some_value = 0;
+
+//[policy_ref_snip1
+
+#include <boost/math/special_functions/gamma.hpp>
+
+using namespace boost::math::policies;
+using namespace boost::math;
+
+// Define a policy:
+typedef policy<
+ domain_error<errno_on_error>,
+ pole_error<errno_on_error>,
+ overflow_error<errno_on_error>,
+ policies::evaluation_error<errno_on_error>
+ > my_policy;
+
+// call the function:
+double t1 = tgamma(some_value, my_policy());
+
+// Alternatively we could use make_policy and define everything at the
call site:
+double t2 = tgamma(some_value, make_policy(
+ domain_error<errno_on_error>(),
+ pole_error<errno_on_error>(),
+ overflow_error<errno_on_error>(),
+ policies::evaluation_error<errno_on_error>()
+ ));
+
+//]
+
+#include <iostream>
+
+int main()
+{
+ std::cout << t1 << " " << t2 << std::endl;
+}
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_ref_snip10.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,25 @@
+// Copyright John Maddock 2007.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Note that this file contains quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+//[policy_ref_snip10
+
+#include <boost/math/special_functions/gamma.hpp>
+
+using namespace boost::math;
+using namespace boost::math::policies;
+
+double t = tgamma(12, policy<digits10<5> >());
+
+//]
+
+#include <iostream>
+
+int main()
+{
+ std::cout << t << std::endl;
+}
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_ref_snip11.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,27 @@
+// Copyright John Maddock 2007.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Note that this file contains quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+//[policy_ref_snip11
+
+#include <boost/math/distributions/normal.hpp>
+
+using namespace boost::math;
+using namespace boost::math::policies;
+
+double q = quantile(
+ normal_distribution<double, policy<digits2<25> > >(),
+ 0.05);
+
+//]
+
+#include <iostream>
+
+int main()
+{
+ std::cout << q << std::endl;
+}
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_ref_snip12.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,38 @@
+// Copyright John Maddock 2007.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Note that this file contains quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+//[policy_ref_snip12
+
+#include <boost/math/special_functions/gamma.hpp>
+
+namespace myspace{
+
+using namespace boost::math::policies;
+
+// Define a policy that does not throw on overflow:
+typedef policy<overflow_error<errno_on_error> > my_policy;
+
+// Define the special functions in this scope to use the policy:
+BOOST_MATH_DECLARE_SPECIAL_FUNCTIONS(my_policy)
+
+}
+
+//
+// Now we can use myspace::tgamma etc.
+// They will automatically use "my_policy":
+//
+double t = myspace::tgamma(30.0); // will not throw on overflow
+
+//]
+
+#include <iostream>
+
+int main()
+{
+ std::cout << t << std::endl;
+}
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_ref_snip13.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,54 @@
+// Copyright John Maddock 2007.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Note that this file contains quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+#ifdef _MSC_VER
+#pragma warning (disable : 4189) // 'd' : local variable is initialized
but not referenced
+#endif
+
+//[policy_ref_snip13
+
+#include <boost/math/distributions/cauchy.hpp>
+
+namespace myspace{
+
+using namespace boost::math::policies;
+
+// Define a policy to use, in this case we want all the distribution
+// accessor functions to compile, even if they are mathematically
+// undefined:
+typedef policy<assert_undefined<false> > my_policy;
+
+BOOST_MATH_DECLARE_DISTRIBUTIONS(double, my_policy)
+
+}
+
+// Now we can use myspace::cauchy etc, which will use policy
+// myspace::mypolicy:
+//
+// This compiles but raises a domain error at runtime:
+//
+void test_cauchy()
+{
+ try
+ {
+ double d = mean(myspace::cauchy());
+ }
+ catch(const std::domain_error& e)
+ {
+ std::cout << e.what() << std::endl;
+ }
+}
+
+//]
+
+#include <iostream>
+
+int main()
+{
+ test_cauchy();
+}
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_ref_snip2.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,34 @@
+// Copyright John Maddock 2007.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Note that this file contains quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+//[policy_ref_snip2
+
+#include <boost/math/distributions/normal.hpp>
+
+using namespace boost::math::policies;
+using namespace boost::math;
+
+// Define a policy:
+typedef policy<
+ overflow_error<ignore_error>
+ > my_policy;
+
+// Define the distribution:
+typedef normal_distribution<double, my_policy> my_norm;
+
+// Get a quantile:
+double q = quantile(my_norm(), 0.05);
+
+//]
+
+#include <iostream>
+
+int main()
+{
+ std::cout << q << std::endl;
+}
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_ref_snip3.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,36 @@
+// Copyright John Maddock 2007.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Note that this file contains quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+double some_value = 2;
+
+//[policy_ref_snip3
+
+#include <boost/math/special_functions/gamma.hpp>
+
+using namespace boost::math::policies;
+using namespace boost::math;
+
+// Define a policy:
+typedef policy<
+ promote_double<false>
+ > my_policy;
+
+// Call the function:
+double t1 = tgamma(some_value, my_policy());
+
+// Alternatively we could use make_policy and define everything at the
call site:
+double t2 = tgamma(some_value, make_policy(promote_double<false>()));
+
+//]
+
+#include <iostream>
+
+int main()
+{
+ std::cout << t1 << " " << t2 << std::endl;
+}
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_ref_snip4.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,40 @@
+// Copyright John Maddock 2007.
+// Copyright Paul A. Bristow 2007.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Note that this file contains quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+#ifdef _MSC_VER
+# pragma warning (disable : 4305) // 'initializing' : truncation
from 'long double' to 'const eval_type'
+# pragma warning (disable : 4244) // 'conversion' : truncation from 'long
double' to 'const eval_type'
+#endif
+
+//[policy_ref_snip4
+
+#include <boost/math/distributions/normal.hpp>
+
+using namespace boost::math::policies;
+using namespace boost::math;
+
+// Define a policy:
+typedef policy<
+ promote_float<false>
+ > my_policy;
+
+// Define the distribution:
+typedef normal_distribution<float, my_policy> my_norm;
+
+// Get a quantile:
+float q = quantile(my_norm(), 0.05f);
+
+//] [policy_ref_snip4]
+
+#include <iostream>
+
+int main()
+{
+ std::cout << q << std::endl;
+}
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_ref_snip5.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,33 @@
+// Copyright John Maddock 2007.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Note that this file contains quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+//[policy_ref_snip5
+
+#include <boost/math/distributions/negative_binomial.hpp>
+
+using namespace boost::math;
+using namespace boost::math::policies;
+
+typedef negative_binomial_distribution<
+ double,
+ policy<discrete_quantile<integer_round_inwards> >
+ > dist_type;
+
+// Lower quantile:
+double x = quantile(dist_type(20, 0.3), 0.05);
+// Upper quantile:
+double y = quantile(complement(dist_type(20, 0.3), 0.05));
+
+//]
+
+#include <iostream>
+
+int main()
+{
+ std::cout << x << " " << y << std::endl;
+}
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_ref_snip6.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,27 @@
+// Copyright John Maddock 2007.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Note that this file contains quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+//[policy_ref_snip6
+
+#include <boost/math/distributions/negative_binomial.hpp>
+
+using namespace boost::math;
+
+// Lower quantile rounded down:
+double x = quantile(negative_binomial(20, 0.3), 0.05);
+// Upper quantile rounded up:
+double y = quantile(complement(negative_binomial(20, 0.3), 0.05));
+
+//]
+
+#include <iostream>
+
+int main()
+{
+ std::cout << x << " " << y << std::endl;
+}
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_ref_snip7.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,33 @@
+// Copyright John Maddock 2007.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Note that this file contains quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+//[policy_ref_snip7
+
+#include <boost/math/distributions/negative_binomial.hpp>
+
+using namespace boost::math;
+using namespace boost::math::policies;
+
+typedef negative_binomial_distribution<
+ double,
+ policy<discrete_quantile<integer_round_inwards> >
+ > dist_type;
+
+// Lower quantile rounded up:
+double x = quantile(dist_type(20, 0.3), 0.05);
+// Upper quantile rounded down:
+double y = quantile(complement(dist_type(20, 0.3), 0.05));
+
+//]
+
+#include <iostream>
+
+int main()
+{
+ std::cout << x << " " << y << std::endl;
+}
=======================================
--- /dev/null
+++ /trunk/libs/math/example/policy_ref_snip8.cpp Wed Jul 29 07:28:12 2009
@@ -0,0 +1,33 @@
+// Copyright John Maddock 2007.
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0. (See accompanying file
+// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Note that this file contains quickbook mark-up as well as code
+// and comments, don't change any of the special comment mark-ups!
+
+//[policy_ref_snip8
+
+#include <boost/math/distributions/negative_binomial.hpp>
+
+using namespace boost::math;
+using namespace boost::math::policies;
+
+typedef negative_binomial_distribution<
+ double,
+ policy<discrete_quantile<integer_round_nearest> >
+ > dist_type;
+
+// Lower quantile rounded up:
+double x = quantile(dist_type(20, 0.3), 0.05);
+// Upper quantile rounded down:
+double y = quantile(complement(dist_type(20, 0.3), 0.05));
+
+//]
+
+#include <iostream>
+
+int main()
+{
+ std::cout << x << " " << y << std::endl;
+}
=======================================
***Additional files exist in this changeset.***
Other related posts:
- » [boost-doc-zh] r283 commited - 校验，补充和修正 - codesite-noreply
