ExRandom  3.0
aux_info.hpp
Go to the documentation of this file.
1 /**
2  * @file aux_info.hpp
3  * @author Charles Karney <charles.karney@sri.com>
4  * @brief Definition of aux_info
5  *
6  * Copyright (c) Charles Karney (2014-2020) and licensed under the MIT/X11
7  * License. For more information, see http://exrandom.sourceforge.net/
8  */
9 
10 #if !defined(EXRANDOM_AUX_INFO_HPP)
11 #define EXRANDOM_AUX_INFO_HPP 1
12 
14 
15 #include <cmath>
16 #include <limits> // For numeric_limits::quiet_NaN
17 
18 namespace exrandom {
19 
20  /**
21  * @brief Some auxiliary functions related to normal and exponential
22  * distributions.
23  *
24  * These are implemented as static functions. These are only needed to
25  * implement the chi-squared tests and to compute the cost and toll of the
26  * algorithms.
27  */
28  class aux_info {
29  public:
30  /**
31  * @return the value of pi.
32  */
33  static double pi() { return std::atan2(0.0, -1.0); }
34 
35  /**
36  * The cumulative uniform function.
37  *
38  * @param x
39  * @return the cumulative probabilty at @e x.
40  */
41  static double cumulative_uniform(double x)
42  { return x < 0 ? 0 : x < 1 ? x : 1; }
43 
44  /**
45  * The cumulative exponential function.
46  *
47  * @param x
48  * @return the cumulative probabilty at @e x.
49  */
50  static double cumulative_exponential(double x) {
51  return x > 0 ? -std::expm1(-x) : 0;
52  }
53 
54  /**
55  * The cumulative normal function.
56  *
57  * @param x
58  * @return the cumulative probabilty at @e x.
59  */
60  static double cumulative_normal(double x)
61  { return std::erf(x / std::sqrt(2.0)) / 2; }
62 
63  /**
64  * @return the entropy of the unit uniform distribution.
65  */
66  static double uniform_entropy() { return 0; }
67 
68  /**
69  * @return the entropy of the unit exponential distribution.
70  */
71  static double exponential_entropy() { return 1; }
72 
73  /**
74  * @return the entropy of the unit exponential distribution.
75  */
76  static double normal_entropy() { return std::log(2 * pi()) / 2 + 0.5; }
77 
78  /**
79  * @return the mean exponent of the floating point representation of a unit
80  * uniform deviate.
81  */
82  static double uniform_mean_exponent() { return -1; }
83 
84  /**
85  * @return the mean exponent of the floating point representation of a unit
86  * exponential deviate.
87  */
88  static double exponential_mean_exponent() {
89  double z = 0;
90  double y0 = cumulative_exponential(0);
91  for (int k = -std::numeric_limits<double>::digits - 4;; ++k) {
92  double y1 = cumulative_exponential(std::pow(2.0, k));
93  z += k * (y1 - y0);
94  if (!(y1 < 1)) break; // Make sure NaN causes an exit
95  y0 = y1;
96  }
97  return z;
98  }
99 
100  /**
101  * @return the mean exponent of the floating point representation of a unit
102  * normal deviate.
103  */
104  static double normal_mean_exponent() {
105  double z = 0;
106  // times 2 to account for negative values too
107  double y0 = 2 * cumulative_normal(0);
108  for (int k = -std::numeric_limits<double>::digits - 4;; ++k) {
109  double y1 = 2 * cumulative_normal(std::pow(2.0, k));
110  z += k * (y1 - y0);
111  if (!(y1 < 1)) break; // Make sure NaN causes an exit
112  y0 = y1;
113  }
114  return z;
115  }
116 
117  /**
118  * @tparam param_type the type of the discrete_normal_dist::param_type.
119  * @param p the parameters for the discrete_normal_dist.
120  * @param i
121  * @param norm the normalizing constant.
122  * @return the probabilty at @e i (unnormalized if norm = 1).
123  */
124  template <typename param_type>
125  static double
126  discrete_normal_prob(const param_type& p, int i, double norm = 1) {
127  double mu = p.mu_num() / double( p.mu_den () ),
128  sigma = p.sigma_num() / double( p.sigma_den () ),
129  x = (i - mu)/sigma;
130  return std::exp( - x * x / 2) / norm;
131  }
132  /**
133  * @tparam param_type the type of the discrete_normal_dist::param_type.
134  * @param p the parameters for the discrete_normal_dist.
135  * @param H the entropy of the distribution.
136  * @return the normalizing constant.
137  */
138  template <typename param_type>
139  static double discrete_normal_norm(const param_type& p, double& H) {
140  int imu = p.mu_num() / p.mu_den(),
141  isig = ( p.sigma_num() + p.sigma_den() - 1 ) / p.sigma_den();
142  double s;
143  if (isig < 10) {
144  s = 0; H = 0;
145  for (int i = imu - 10 * isig; i < imu + 10 * isig; ++i) {
146  double P = discrete_normal_prob(p, i, 1.0);
147  s += P;
148  if (P > 0)
149  H -= P * std::log(P);
150  }
151  H = H/s + std::log(s);
152  } else {
153  // Continuous approximations are very good
154  double sigma = p.sigma_num() / double(p.sigma_den());
155  s = std::sqrt(2 * pi()) * sigma;
156  H = std::log(s) + 0.5;
157  }
158  return s;
159  }
160  };
161 
162 }
163 
164 #endif // EXRANDOM_AUX_INFO_HPP
static double uniform_mean_exponent()
Definition: aux_info.hpp:82
static double discrete_normal_norm(const param_type &p, double &H)
Definition: aux_info.hpp:139
static double uniform_entropy()
Definition: aux_info.hpp:66
static double cumulative_normal(double x)
Definition: aux_info.hpp:60
static double cumulative_exponential(double x)
Definition: aux_info.hpp:50
static double normal_mean_exponent()
Definition: aux_info.hpp:104
static double cumulative_uniform(double x)
Definition: aux_info.hpp:41
Configuration for exrandom.
static double exponential_mean_exponent()
Definition: aux_info.hpp:88
static double pi()
Definition: aux_info.hpp:33
static double discrete_normal_prob(const param_type &p, int i, double norm=1)
Definition: aux_info.hpp:126
The common namespace.
Definition: aux_info.hpp:18
Some auxiliary functions related to normal and exponential distributions.
Definition: aux_info.hpp:28
static double normal_entropy()
Definition: aux_info.hpp:76
static double exponential_entropy()
Definition: aux_info.hpp:71