To check out this repository please hg clone the following URL, or open the URL using EasyMercurial or your preferred Mercurial client.

The primary repository for this project is hosted at https://github.com/sonic-visualiser/sv-dependency-builds .
This repository is a read-only copy which is updated automatically every hour.

Statistics Download as Zip
| Branch: | Tag: | Revision:

root / any / include / boost / math / distributions / triangular.hpp @ 160:cff480c41f97

History | View | Annotate | Download (18.3 KB)

1
//  Copyright John Maddock 2006, 2007.
2
//  Copyright Paul A. Bristow 2006, 2007.
3
//  Use, modification and distribution are subject to the
4
//  Boost Software License, Version 1.0. (See accompanying file
5
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
6

    
7
#ifndef BOOST_STATS_TRIANGULAR_HPP
8
#define BOOST_STATS_TRIANGULAR_HPP
9

    
10
// http://mathworld.wolfram.com/TriangularDistribution.html
11
// Note that the 'constructors' defined by Wolfram are difference from those here,
12
// for example
13
// N[variance[triangulardistribution{1, +2}, 1.5], 50] computes 
14
// 0.041666666666666666666666666666666666666666666666667
15
// TriangularDistribution{1, +2}, 1.5 is the analog of triangular_distribution(1, 1.5, 2)
16

    
17
// http://en.wikipedia.org/wiki/Triangular_distribution
18

    
19
#include <boost/math/distributions/fwd.hpp>
20
#include <boost/math/special_functions/expm1.hpp>
21
#include <boost/math/distributions/detail/common_error_handling.hpp>
22
#include <boost/math/distributions/complement.hpp>
23
#include <boost/math/constants/constants.hpp>
24

    
25
#include <utility>
26

    
27
namespace boost{ namespace math
28
{
29
  namespace detail
30
  {
31
    template <class RealType, class Policy>
32
    inline bool check_triangular_lower(
33
      const char* function,
34
      RealType lower,
35
      RealType* result, const Policy& pol)
36
    {
37
      if((boost::math::isfinite)(lower))
38
      { // Any finite value is OK.
39
        return true;
40
      }
41
      else
42
      { // Not finite: infinity or NaN.
43
        *result = policies::raise_domain_error<RealType>(
44
          function,
45
          "Lower parameter is %1%, but must be finite!", lower, pol);
46
        return false;
47
      }
48
    } // bool check_triangular_lower(
49

    
50
    template <class RealType, class Policy>
51
    inline bool check_triangular_mode(
52
      const char* function,
53
      RealType mode,
54
      RealType* result, const Policy& pol)
55
    {
56
      if((boost::math::isfinite)(mode))
57
      { // any finite value is OK.
58
        return true;
59
      }
60
      else
61
      { // Not finite: infinity or NaN.
62
        *result = policies::raise_domain_error<RealType>(
63
          function,
64
          "Mode parameter is %1%, but must be finite!", mode, pol);
65
        return false;
66
      }
67
    } // bool check_triangular_mode(
68

    
69
    template <class RealType, class Policy>
70
    inline bool check_triangular_upper(
71
      const char* function,
72
      RealType upper,
73
      RealType* result, const Policy& pol)
74
    {
75
      if((boost::math::isfinite)(upper))
76
      { // any finite value is OK.
77
        return true;
78
      }
79
      else
80
      { // Not finite: infinity or NaN.
81
        *result = policies::raise_domain_error<RealType>(
82
          function,
83
          "Upper parameter is %1%, but must be finite!", upper, pol);
84
        return false;
85
      }
86
    } // bool check_triangular_upper(
87

    
88
    template <class RealType, class Policy>
89
    inline bool check_triangular_x(
90
      const char* function,
91
      RealType const& x,
92
      RealType* result, const Policy& pol)
93
    {
94
      if((boost::math::isfinite)(x))
95
      { // Any finite value is OK
96
        return true;
97
      }
98
      else
99
      { // Not finite: infinity or NaN.
100
        *result = policies::raise_domain_error<RealType>(
101
          function,
102
          "x parameter is %1%, but must be finite!", x, pol);
103
        return false;
104
      }
105
    } // bool check_triangular_x
106

    
107
    template <class RealType, class Policy>
108
    inline bool check_triangular(
109
      const char* function,
110
      RealType lower,
111
      RealType mode,
112
      RealType upper,
113
      RealType* result, const Policy& pol)
114
    {
115
      if ((check_triangular_lower(function, lower, result, pol) == false)
116
        || (check_triangular_mode(function, mode, result, pol) == false)
117
        || (check_triangular_upper(function, upper, result, pol) == false))
118
      { // Some parameter not finite.
119
        return false;
120
      }
121
      else if (lower >= upper) // lower == upper NOT useful.
122
      { // lower >= upper.
123
        *result = policies::raise_domain_error<RealType>(
124
          function,
125
          "lower parameter is %1%, but must be less than upper!", lower, pol);
126
        return false;
127
      }
128
      else
129
      { // Check lower <= mode <= upper.
130
        if (mode < lower)
131
        {
132
          *result = policies::raise_domain_error<RealType>(
133
            function,
134
            "mode parameter is %1%, but must be >= than lower!", lower, pol);
135
          return false;
136
        }
137
        if (mode > upper)
138
        {
139
          *result = policies::raise_domain_error<RealType>(
140
            function,
141
            "mode parameter is %1%, but must be <= than upper!", upper, pol);
142
          return false;
143
        }
144
        return true; // All OK.
145
      }
146
    } // bool check_triangular
147
  } // namespace detail
148

    
149
  template <class RealType = double, class Policy = policies::policy<> >
150
  class triangular_distribution
151
  {
152
  public:
153
    typedef RealType value_type;
154
    typedef Policy policy_type;
155

    
156
    triangular_distribution(RealType l_lower = -1, RealType l_mode = 0, RealType l_upper = 1)
157
      : m_lower(l_lower), m_mode(l_mode), m_upper(l_upper) // Constructor.
158
    { // Evans says 'standard triangular' is lower 0, mode 1/2, upper 1,
159
      // has median sqrt(c/2) for c <=1/2 and 1 - sqrt(1-c)/2 for c >= 1/2
160
      // But this -1, 0, 1 is more useful in most applications to approximate normal distribution,
161
      // where the central value is the most likely and deviations either side equally likely.
162
      RealType result;
163
      detail::check_triangular("boost::math::triangular_distribution<%1%>::triangular_distribution",l_lower, l_mode, l_upper, &result, Policy());
164
    }
165
    // Accessor functions.
166
    RealType lower()const
167
    {
168
      return m_lower;
169
    }
170
    RealType mode()const
171
    {
172
      return m_mode;
173
    }
174
    RealType upper()const
175
    {
176
      return m_upper;
177
    }
178
  private:
179
    // Data members:
180
    RealType m_lower;  // distribution lower aka a
181
    RealType m_mode;  // distribution mode aka c
182
    RealType m_upper;  // distribution upper aka b
183
  }; // class triangular_distribution
184

    
185
  typedef triangular_distribution<double> triangular;
186

    
187
  template <class RealType, class Policy>
188
  inline const std::pair<RealType, RealType> range(const triangular_distribution<RealType, Policy>& /* dist */)
189
  { // Range of permissible values for random variable x.
190
    using boost::math::tools::max_value;
191
    return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>());
192
  }
193

    
194
  template <class RealType, class Policy>
195
  inline const std::pair<RealType, RealType> support(const triangular_distribution<RealType, Policy>& dist)
196
  { // Range of supported values for random variable x.
197
    // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
198
    return std::pair<RealType, RealType>(dist.lower(), dist.upper());
199
  }
200

    
201
  template <class RealType, class Policy>
202
  RealType pdf(const triangular_distribution<RealType, Policy>& dist, const RealType& x)
203
  {
204
    static const char* function = "boost::math::pdf(const triangular_distribution<%1%>&, %1%)";
205
    RealType lower = dist.lower();
206
    RealType mode = dist.mode();
207
    RealType upper = dist.upper();
208
    RealType result = 0; // of checks.
209
    if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy()))
210
    {
211
      return result;
212
    }
213
    if(false == detail::check_triangular_x(function, x, &result, Policy()))
214
    {
215
      return result;
216
    }
217
    if((x < lower) || (x > upper))
218
    {
219
      return 0;
220
    }
221
    if (x == lower)
222
    { // (mode - lower) == 0 which would lead to divide by zero!
223
      return (mode == lower) ? 2 / (upper - lower) : RealType(0);
224
    }
225
    else if (x == upper)
226
    {
227
      return (mode == upper) ? 2 / (upper - lower) : RealType(0);
228
    }
229
    else if (x <= mode)
230
    {
231
      return 2 * (x - lower) / ((upper - lower) * (mode - lower));
232
    }
233
    else
234
    {  // (x > mode)
235
      return 2 * (upper - x) / ((upper - lower) * (upper - mode));
236
    }
237
  } // RealType pdf(const triangular_distribution<RealType, Policy>& dist, const RealType& x)
238

    
239
  template <class RealType, class Policy>
240
  inline RealType cdf(const triangular_distribution<RealType, Policy>& dist, const RealType& x)
241
  {
242
    static const char* function = "boost::math::cdf(const triangular_distribution<%1%>&, %1%)";
243
    RealType lower = dist.lower();
244
    RealType mode = dist.mode();
245
    RealType upper = dist.upper();
246
    RealType result = 0; // of checks.
247
    if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy()))
248
    {
249
      return result;
250
    }
251
    if(false == detail::check_triangular_x(function, x, &result, Policy()))
252
    {
253
      return result;
254
    }
255
    if((x <= lower))
256
    {
257
      return 0;
258
    }
259
    if (x >= upper)
260
    {
261
      return 1;
262
    }
263
    // else lower < x < upper
264
    if (x <= mode)
265
    {
266
      return ((x - lower) * (x - lower)) / ((upper - lower) * (mode - lower));
267
    }
268
    else
269
    {
270
      return 1 - (upper - x) *  (upper - x) / ((upper - lower) * (upper - mode));
271
    }
272
  } // RealType cdf(const triangular_distribution<RealType, Policy>& dist, const RealType& x)
273

    
274
  template <class RealType, class Policy>
275
  RealType quantile(const triangular_distribution<RealType, Policy>& dist, const RealType& p)
276
  {
277
    BOOST_MATH_STD_USING  // for ADL of std functions (sqrt).
278
    static const char* function = "boost::math::quantile(const triangular_distribution<%1%>&, %1%)";
279
    RealType lower = dist.lower();
280
    RealType mode = dist.mode();
281
    RealType upper = dist.upper();
282
    RealType result = 0; // of checks
283
    if(false == detail::check_triangular(function,lower, mode, upper, &result, Policy()))
284
    {
285
      return result;
286
    }
287
    if(false == detail::check_probability(function, p, &result, Policy()))
288
    {
289
      return result;
290
    }
291
    if(p == 0)
292
    {
293
      return lower;
294
    }
295
    if(p == 1)
296
    {
297
      return upper;
298
    }
299
    RealType p0 = (mode - lower) / (upper - lower);
300
    RealType q = 1 - p;
301
    if (p < p0)
302
    {
303
      result = sqrt((upper - lower) * (mode - lower) * p) + lower;
304
    }
305
    else if (p == p0)
306
    {
307
      result = mode;
308
    }
309
    else // p > p0
310
    {
311
      result = upper - sqrt((upper - lower) * (upper - mode) * q);
312
    }
313
    return result;
314

    
315
  } // RealType quantile(const triangular_distribution<RealType, Policy>& dist, const RealType& q)
316

    
317
  template <class RealType, class Policy>
318
  RealType cdf(const complemented2_type<triangular_distribution<RealType, Policy>, RealType>& c)
319
  {
320
    static const char* function = "boost::math::cdf(const triangular_distribution<%1%>&, %1%)";
321
    RealType lower = c.dist.lower();
322
    RealType mode = c.dist.mode();
323
    RealType upper = c.dist.upper();
324
    RealType x = c.param;
325
    RealType result = 0; // of checks.
326
    if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy()))
327
    {
328
      return result;
329
    }
330
    if(false == detail::check_triangular_x(function, x, &result, Policy()))
331
    {
332
      return result;
333
    }
334
    if (x <= lower)
335
    {
336
      return 1;
337
    }
338
    if (x >= upper)
339
    {
340
      return 0;
341
    }
342
    if (x <= mode)
343
    {
344
      return 1 - ((x - lower) * (x - lower)) / ((upper - lower) * (mode - lower));
345
    }
346
    else
347
    {
348
      return (upper - x) *  (upper - x) / ((upper - lower) * (upper - mode));
349
    }
350
  } // RealType cdf(const complemented2_type<triangular_distribution<RealType, Policy>, RealType>& c)
351

    
352
  template <class RealType, class Policy>
353
  RealType quantile(const complemented2_type<triangular_distribution<RealType, Policy>, RealType>& c)
354
  {
355
    BOOST_MATH_STD_USING  // Aid ADL for sqrt.
356
    static const char* function = "boost::math::quantile(const triangular_distribution<%1%>&, %1%)";
357
    RealType l = c.dist.lower();
358
    RealType m = c.dist.mode();
359
    RealType u = c.dist.upper();
360
    RealType q = c.param; // probability 0 to 1.
361
    RealType result = 0; // of checks.
362
    if(false == detail::check_triangular(function, l, m, u, &result, Policy()))
363
    {
364
      return result;
365
    }
366
    if(false == detail::check_probability(function, q, &result, Policy()))
367
    {
368
      return result;
369
    }
370
    if(q == 0)
371
    {
372
      return u;
373
    }
374
    if(q == 1)
375
    {
376
      return l;
377
    }
378
    RealType lower = c.dist.lower();
379
    RealType mode = c.dist.mode();
380
    RealType upper = c.dist.upper();
381

    
382
    RealType p = 1 - q;
383
    RealType p0 = (mode - lower) / (upper - lower);
384
    if(p < p0)
385
    {
386
      RealType s = (upper - lower) * (mode - lower);
387
      s *= p;
388
      result = sqrt((upper - lower) * (mode - lower) * p) + lower;
389
    }
390
    else if (p == p0)
391
    {
392
      result = mode;
393
    }
394
    else // p > p0
395
    {
396
      result = upper - sqrt((upper - lower) * (upper - mode) * q);
397
    }
398
    return result;
399
  } // RealType quantile(const complemented2_type<triangular_distribution<RealType, Policy>, RealType>& c)
400

    
401
  template <class RealType, class Policy>
402
  inline RealType mean(const triangular_distribution<RealType, Policy>& dist)
403
  {
404
    static const char* function = "boost::math::mean(const triangular_distribution<%1%>&)";
405
    RealType lower = dist.lower();
406
    RealType mode = dist.mode();
407
    RealType upper = dist.upper();
408
    RealType result = 0;  // of checks.
409
    if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy()))
410
    {
411
      return result;
412
    }
413
    return (lower + upper + mode) / 3;
414
  } // RealType mean(const triangular_distribution<RealType, Policy>& dist)
415

    
416

    
417
  template <class RealType, class Policy>
418
  inline RealType variance(const triangular_distribution<RealType, Policy>& dist)
419
  {
420
    static const char* function = "boost::math::mean(const triangular_distribution<%1%>&)";
421
    RealType lower = dist.lower();
422
    RealType mode = dist.mode();
423
    RealType upper = dist.upper();
424
    RealType result = 0; // of checks.
425
    if(false == detail::check_triangular(function, lower, mode, upper, &result, Policy()))
426
    {
427
      return result;
428
    }
429
    return (lower * lower + upper * upper + mode * mode - lower * upper - lower * mode - upper * mode) / 18;
430
  } // RealType variance(const triangular_distribution<RealType, Policy>& dist)
431

    
432
  template <class RealType, class Policy>
433
  inline RealType mode(const triangular_distribution<RealType, Policy>& dist)
434
  {
435
    static const char* function = "boost::math::mode(const triangular_distribution<%1%>&)";
436
    RealType mode = dist.mode();
437
    RealType result = 0; // of checks.
438
    if(false == detail::check_triangular_mode(function, mode, &result, Policy()))
439
    { // This should never happen!
440
      return result;
441
    }
442
    return mode;
443
  } // RealType mode
444

    
445
  template <class RealType, class Policy>
446
  inline RealType median(const triangular_distribution<RealType, Policy>& dist)
447
  {
448
    BOOST_MATH_STD_USING // ADL of std functions.
449
    static const char* function = "boost::math::median(const triangular_distribution<%1%>&)";
450
    RealType mode = dist.mode();
451
    RealType result = 0; // of checks.
452
    if(false == detail::check_triangular_mode(function, mode, &result, Policy()))
453
    { // This should never happen!
454
      return result;
455
    }
456
    RealType lower = dist.lower();
457
    RealType upper = dist.upper();
458
    if (mode >= (upper + lower) / 2)
459
    {
460
      return lower + sqrt((upper - lower) * (mode - lower)) / constants::root_two<RealType>();
461
    }
462
    else
463
    {
464
      return upper - sqrt((upper - lower) * (upper - mode)) / constants::root_two<RealType>();
465
    }
466
  } // RealType mode
467

    
468
  template <class RealType, class Policy>
469
  inline RealType skewness(const triangular_distribution<RealType, Policy>& dist)
470
  {
471
    BOOST_MATH_STD_USING  // for ADL of std functions
472
    using namespace boost::math::constants; // for root_two
473
    static const char* function = "boost::math::skewness(const triangular_distribution<%1%>&)";
474

    
475
    RealType lower = dist.lower();
476
    RealType mode = dist.mode();
477
    RealType upper = dist.upper();
478
    RealType result = 0; // of checks.
479
    if(false == boost::math::detail::check_triangular(function,lower, mode, upper, &result, Policy()))
480
    {
481
      return result;
482
    }
483
    return root_two<RealType>() * (lower + upper - 2 * mode) * (2 * lower - upper - mode) * (lower - 2 * upper + mode) /
484
      (5 * pow((lower * lower + upper * upper + mode * mode 
485
        - lower * upper - lower * mode - upper * mode), RealType(3)/RealType(2)));
486
    // #11768: Skewness formula for triangular distribution is incorrect -  corrected 29 Oct 2015 for release 1.61.
487
  } // RealType skewness(const triangular_distribution<RealType, Policy>& dist)
488

    
489
  template <class RealType, class Policy>
490
  inline RealType kurtosis(const triangular_distribution<RealType, Policy>& dist)
491
  { // These checks may be belt and braces as should have been checked on construction?
492
    static const char* function = "boost::math::kurtosis(const triangular_distribution<%1%>&)";
493
    RealType lower = dist.lower();
494
    RealType upper = dist.upper();
495
    RealType mode = dist.mode();
496
    RealType result = 0;  // of checks.
497
    if(false == detail::check_triangular(function,lower, mode, upper, &result, Policy()))
498
    {
499
      return result;
500
    }
501
    return static_cast<RealType>(12)/5; //  12/5 = 2.4;
502
  } // RealType kurtosis_excess(const triangular_distribution<RealType, Policy>& dist)
503

    
504
  template <class RealType, class Policy>
505
  inline RealType kurtosis_excess(const triangular_distribution<RealType, Policy>& dist)
506
  { // These checks may be belt and braces as should have been checked on construction?
507
    static const char* function = "boost::math::kurtosis_excess(const triangular_distribution<%1%>&)";
508
    RealType lower = dist.lower();
509
    RealType upper = dist.upper();
510
    RealType mode = dist.mode();
511
    RealType result = 0;  // of checks.
512
    if(false == detail::check_triangular(function,lower, mode, upper, &result, Policy()))
513
    {
514
      return result;
515
    }
516
    return static_cast<RealType>(-3)/5; // - 3/5 = -0.6
517
    // Assuming mathworld really means kurtosis excess?  Wikipedia now corrected to match this.
518
  }
519

    
520
} // namespace math
521
} // namespace boost
522

    
523
// This include must be at the end, *after* the accessors
524
// for this distribution have been defined, in order to
525
// keep compilers that support two-phase lookup happy.
526
#include <boost/math/distributions/detail/derived_accessors.hpp>
527

    
528
#endif // BOOST_STATS_TRIANGULAR_HPP
529

    
530

    
531