~tsarev/boostdc/cmake

« back to all changes in this revision

Viewing changes to boost/boost/math/distributions/non_central_chi_squared.hpp

  • Committer: bigmuscle
  • Date: 2010-05-08 08:47:15 UTC
  • Revision ID: svn-v4:5fb55d53-692c-0410-a46a-e90ab66e00ee:trunk:497
removed old boost version

Show diffs side-by-side

added added

removed removed

Lines of Context:
1
 
// boost\math\distributions\non_central_chi_squared.hpp
2
 
 
3
 
// Copyright John Maddock 2008.
4
 
 
5
 
// Use, modification and distribution are subject to the
6
 
// Boost Software License, Version 1.0.
7
 
// (See accompanying file LICENSE_1_0.txt
8
 
// or copy at http://www.boost.org/LICENSE_1_0.txt)
9
 
 
10
 
#ifndef BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP
11
 
#define BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP
12
 
 
13
 
#include <boost/math/distributions/fwd.hpp>
14
 
#include <boost/math/special_functions/gamma.hpp> // for incomplete gamma. gamma_q
15
 
#include <boost/math/special_functions/bessel.hpp> // for cyl_bessel_i
16
 
#include <boost/math/special_functions/round.hpp> // for iround
17
 
#include <boost/math/distributions/complement.hpp> // complements
18
 
#include <boost/math/distributions/chi_squared.hpp> // central distribution
19
 
#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
20
 
#include <boost/math/special_functions/fpclassify.hpp> // isnan.
21
 
#include <boost/math/tools/roots.hpp> // for root finding.
22
 
#include <boost/math/distributions/detail/generic_mode.hpp>
23
 
#include <boost/math/distributions/detail/generic_quantile.hpp>
24
 
 
25
 
namespace boost
26
 
{
27
 
   namespace math
28
 
   {
29
 
 
30
 
      template <class RealType, class Policy>
31
 
      class non_central_chi_squared_distribution;
32
 
 
33
 
      namespace detail{
34
 
 
35
 
         template <class T, class Policy>
36
 
         T non_central_chi_square_q(T x, T f, T theta, const Policy& pol, T init_sum = 0)
37
 
         {
38
 
            //
39
 
            // Computes the complement of the Non-Central Chi-Square
40
 
            // Distribution CDF by summing a weighted sum of complements
41
 
            // of the central-distributions.  The weighting factor is
42
 
            // a Poisson Distribution.
43
 
            //
44
 
            // This is an application of the technique described in:
45
 
            //
46
 
            // Computing discrete mixtures of continuous
47
 
            // distributions: noncentral chisquare, noncentral t
48
 
            // and the distribution of the square of the sample
49
 
            // multiple correlation coeficient.
50
 
            // D. Benton, K. Krishnamoorthy.
51
 
            // Computational Statistics & Data Analysis 43 (2003) 249 - 267
52
 
            //
53
 
            BOOST_MATH_STD_USING
54
 
 
55
 
            // Special case:
56
 
            if(x == 0)
57
 
               return 1;
58
 
 
59
 
            //
60
 
            // Initialize the variables we'll be using:
61
 
            //
62
 
            T lambda = theta / 2;
63
 
            T del = f / 2;
64
 
            T y = x / 2;
65
 
            boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
66
 
            T errtol = boost::math::policies::get_epsilon<T, Policy>();
67
 
            T sum = init_sum;
68
 
            //
69
 
            // k is the starting location for iteration, we'll
70
 
            // move both forwards and backwards from this point.
71
 
            // k is chosen as the peek of the Poisson weights, which
72
 
            // will occur *before* the largest term.
73
 
            //
74
 
            int k = iround(lambda, pol);
75
 
            // Forwards and backwards Poisson weights:
76
 
            T poisf = boost::math::gamma_p_derivative(1 + k, lambda, pol);
77
 
            T poisb = poisf * k / lambda;
78
 
            // Initial forwards central chi squared term:
79
 
            T gamf = boost::math::gamma_q(del + k, y, pol);
80
 
            // Forwards and backwards recursion terms on the central chi squared:
81
 
            T xtermf = boost::math::gamma_p_derivative(del + 1 + k, y, pol);
82
 
            T xtermb = xtermf * (del + k) / y;
83
 
            // Initial backwards central chi squared term:
84
 
            T gamb = gamf - xtermb;
85
 
 
86
 
            //
87
 
            // Forwards iteration first, this is the
88
 
            // stable direction for the gamma function
89
 
            // recurrences:
90
 
            //
91
 
            int i;
92
 
            for(i = k; static_cast<boost::uintmax_t>(i-k) < max_iter; ++i)
93
 
            {
94
 
               T term = poisf * gamf;
95
 
               sum += term;
96
 
               poisf *= lambda / (i + 1);
97
 
               gamf += xtermf;
98
 
               xtermf *= y / (del + i + 1);
99
 
               if(((sum == 0) || (fabs(term / sum) < errtol)) && (term >= poisf * gamf))
100
 
                  break;
101
 
            }
102
 
            //Error check:
103
 
            if(static_cast<boost::uintmax_t>(i-k) >= max_iter)
104
 
               policies::raise_evaluation_error(
105
 
                  "cdf(non_central_chi_squared_distribution<%1%>, %1%)", 
106
 
                  "Series did not converge, closest value was %1%", sum, pol);
107
 
            //
108
 
            // Now backwards iteration: the gamma
109
 
            // function recurrences are unstable in this
110
 
            // direction, we rely on the terms deminishing in size
111
 
            // faster than we introduce cancellation errors.  
112
 
            // For this reason it's very important that we start
113
 
            // *before* the largest term so that backwards iteration
114
 
            // is strictly converging.
115
 
            //
116
 
            for(i = k - 1; i >= 0; --i)
117
 
            {
118
 
               T term = poisb * gamb;
119
 
               sum += term;
120
 
               poisb *= i / lambda;
121
 
               xtermb *= (del + i) / y;
122
 
               gamb -= xtermb;
123
 
               if((sum == 0) || (fabs(term / sum) < errtol))
124
 
                  break;
125
 
            }
126
 
 
127
 
            return sum;
128
 
         }
129
 
 
130
 
         template <class T, class Policy>
131
 
         T non_central_chi_square_p_ding(T x, T f, T theta, const Policy& pol, T init_sum = 0)
132
 
         {
133
 
            //
134
 
            // This is an implementation of:
135
 
            //
136
 
            // Algorithm AS 275: 
137
 
            // Computing the Non-Central #2 Distribution Function
138
 
            // Cherng G. Ding
139
 
            // Applied Statistics, Vol. 41, No. 2. (1992), pp. 478-482.
140
 
            //
141
 
            // This uses a stable forward iteration to sum the
142
 
            // CDF, unfortunately this can not be used for large
143
 
            // values of the non-centrality parameter because:
144
 
            // * The first term may underfow to zero.
145
 
            // * We may need an extra-ordinary number of terms
146
 
            //   before we reach the first *significant* term.
147
 
            //
148
 
            BOOST_MATH_STD_USING
149
 
            // Special case:
150
 
            if(x == 0)
151
 
               return 0;
152
 
            T tk = boost::math::gamma_p_derivative(f/2 + 1, x/2, pol);
153
 
            T lambda = theta / 2;
154
 
            T vk = exp(-lambda);
155
 
            T uk = vk;
156
 
            T sum = init_sum + tk * vk;
157
 
            if(sum == 0)
158
 
               return sum;
159
 
 
160
 
            boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
161
 
            T errtol = boost::math::policies::get_epsilon<T, Policy>();
162
 
 
163
 
            int i;
164
 
            T lterm(0), term(0);
165
 
            for(i = 1; static_cast<boost::uintmax_t>(i) < max_iter; ++i)
166
 
            {
167
 
               tk = tk * x / (f + 2 * i);
168
 
               uk = uk * lambda / i;
169
 
               vk = vk + uk;
170
 
               lterm = term;
171
 
               term = vk * tk;
172
 
               sum += term;
173
 
               if((fabs(term / sum) < errtol) && (term <= lterm))
174
 
                  break;
175
 
            }
176
 
            //Error check:
177
 
            if(static_cast<boost::uintmax_t>(i) >= max_iter)
178
 
               policies::raise_evaluation_error(
179
 
                  "cdf(non_central_chi_squared_distribution<%1%>, %1%)", 
180
 
                  "Series did not converge, closest value was %1%", sum, pol);
181
 
            return sum;
182
 
         }
183
 
 
184
 
 
185
 
         template <class T, class Policy>
186
 
         T non_central_chi_square_p(T y, T n, T lambda, const Policy& pol, T init_sum)
187
 
         {
188
 
            //
189
 
            // This is taken more or less directly from:
190
 
            //
191
 
            // Computing discrete mixtures of continuous
192
 
            // distributions: noncentral chisquare, noncentral t
193
 
            // and the distribution of the square of the sample
194
 
            // multiple correlation coeficient.
195
 
            // D. Benton, K. Krishnamoorthy.
196
 
            // Computational Statistics & Data Analysis 43 (2003) 249 - 267
197
 
            //
198
 
            // We're summing a Poisson weighting term multiplied by
199
 
            // a central chi squared distribution.
200
 
            //
201
 
            BOOST_MATH_STD_USING
202
 
            // Special case:
203
 
            if(y == 0)
204
 
               return 0;
205
 
            boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
206
 
            T errtol = boost::math::policies::get_epsilon<T, Policy>();
207
 
            T errorf(0), errorb(0);
208
 
 
209
 
            T x = y / 2;
210
 
            T del = lambda / 2;
211
 
            //
212
 
            // Starting location for the iteration, we'll iterate
213
 
            // both forwards and backwards from this point.  The
214
 
            // location chosen is the maximum of the Poisson weight
215
 
            // function, which ocurrs *after* the largest term in the
216
 
            // sum.
217
 
            //
218
 
            int k = iround(del, pol);
219
 
            T a = n / 2 + k;
220
 
            // Central chi squared term for forward iteration:
221
 
            T gamkf = boost::math::gamma_p(a, x, pol);
222
 
 
223
 
            if(lambda == 0)
224
 
               return gamkf;
225
 
            // Central chi squared term for backward iteration:
226
 
            T gamkb = gamkf;
227
 
            // Forwards Poisson weight:
228
 
            T poiskf = gamma_p_derivative(k+1, del, pol);
229
 
            // Backwards Poisson weight:
230
 
            T poiskb = poiskf;
231
 
            // Forwards gamma function recursion term:
232
 
            T xtermf = boost::math::gamma_p_derivative(a, x, pol);
233
 
            // Backwards gamma function recursion term:
234
 
            T xtermb = xtermf * x / a;
235
 
            T sum = init_sum + poiskf * gamkf;
236
 
            if(sum == 0)
237
 
               return sum;
238
 
            int i = 1;
239
 
            //
240
 
            // Backwards recursion first, this is the stable
241
 
            // direction for gamma function recurrences:
242
 
            //
243
 
            while(i <= k) 
244
 
            {
245
 
               xtermb *= (a - i + 1) / x;
246
 
               gamkb += xtermb;
247
 
               poiskb = poiskb * (k - i + 1) / del;
248
 
               errorf = errorb;
249
 
               errorb = gamkb * poiskb;
250
 
               sum += errorb;
251
 
               if((fabs(errorb / sum) < errtol) && (errorb <= errorf))
252
 
                  break;
253
 
               ++i;
254
 
            }
255
 
            i = 1;
256
 
            //
257
 
            // Now forwards recursion, the gamma function
258
 
            // recurrence relation is unstable in this direction,
259
 
            // so we rely on the magnitude of successive terms
260
 
            // decreasing faster than we introduce cancellation error.
261
 
            // For this reason it's vital that k is chosen to be *after*
262
 
            // the largest term, so that successive forward iterations
263
 
            // are strictly (and rapidly) converging.
264
 
            //
265
 
            do
266
 
            {
267
 
               xtermf = xtermf * x / (a + i - 1);
268
 
               gamkf = gamkf - xtermf;
269
 
               poiskf = poiskf * del / (k + i);
270
 
               errorf = poiskf * gamkf;
271
 
               sum += errorf;
272
 
               ++i;
273
 
            }while((fabs(errorf / sum) > errtol) && (static_cast<boost::uintmax_t>(i) < max_iter));
274
 
 
275
 
            //Error check:
276
 
            if(static_cast<boost::uintmax_t>(i) >= max_iter)
277
 
               policies::raise_evaluation_error(
278
 
                  "cdf(non_central_chi_squared_distribution<%1%>, %1%)", 
279
 
                  "Series did not converge, closest value was %1%", sum, pol);
280
 
 
281
 
            return sum;
282
 
         }
283
 
 
284
 
         template <class T, class Policy>
285
 
         T non_central_chi_square_pdf(T x, T n, T lambda, const Policy& pol)
286
 
         {
287
 
            //
288
 
            // As above but for the PDF:
289
 
            //
290
 
            BOOST_MATH_STD_USING
291
 
            boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
292
 
            T errtol = boost::math::policies::get_epsilon<T, Policy>();
293
 
            T x2 = x / 2;
294
 
            T n2 = n / 2;
295
 
            T l2 = lambda / 2;
296
 
            T sum = 0;
297
 
            int k = itrunc(l2);
298
 
            T pois = gamma_p_derivative(k + 1, l2, pol) * gamma_p_derivative(n2 + k, x2);
299
 
            if(pois == 0)
300
 
               return 0;
301
 
            T poisb = pois;
302
 
            for(int i = k; ; ++i)
303
 
            {
304
 
               sum += pois;
305
 
               if(pois / sum < errtol)
306
 
                  break;
307
 
               if(static_cast<boost::uintmax_t>(i - k) >= max_iter)
308
 
                  return policies::raise_evaluation_error(
309
 
                     "pdf(non_central_chi_squared_distribution<%1%>, %1%)", 
310
 
                     "Series did not converge, closest value was %1%", sum, pol);
311
 
               pois *= l2 * x2 / ((i + 1) * (n2 + i));
312
 
            }
313
 
            for(int i = k - 1; i >= 0; --i)
314
 
            {
315
 
               poisb *= (i + 1) * (n2 + i) / (l2 * x2);
316
 
               sum += poisb;
317
 
               if(poisb / sum < errtol)
318
 
                  break;
319
 
            }
320
 
            return sum / 2;
321
 
         }
322
 
 
323
 
         template <class RealType, class Policy>
324
 
         inline RealType non_central_chi_squared_cdf(RealType x, RealType k, RealType l, bool invert, const Policy&)
325
 
         {
326
 
            typedef typename policies::evaluation<RealType, Policy>::type value_type;
327
 
            typedef typename policies::normalise<
328
 
               Policy, 
329
 
               policies::promote_float<false>, 
330
 
               policies::promote_double<false>, 
331
 
               policies::discrete_quantile<>,
332
 
               policies::assert_undefined<> >::type forwarding_policy;
333
 
 
334
 
            BOOST_MATH_STD_USING
335
 
            value_type result;
336
 
            if(l == 0)
337
 
               result = cdf(boost::math::chi_squared_distribution<RealType, Policy>(k), x);
338
 
            else if(x > k + l)
339
 
            {
340
 
               // Complement is the smaller of the two:
341
 
               result = detail::non_central_chi_square_q(
342
 
                  static_cast<value_type>(x), 
343
 
                  static_cast<value_type>(k), 
344
 
                  static_cast<value_type>(l), 
345
 
                  forwarding_policy(), 
346
 
                  static_cast<value_type>(invert ? 0 : -1));
347
 
               invert = !invert;
348
 
            }
349
 
            else if(l < 200)
350
 
            {
351
 
               // For small values of the non-centrality parameter
352
 
               // we can use Ding's method:
353
 
               result = detail::non_central_chi_square_p_ding(
354
 
                  static_cast<value_type>(x), 
355
 
                  static_cast<value_type>(k), 
356
 
                  static_cast<value_type>(l), 
357
 
                  forwarding_policy(),
358
 
                  static_cast<value_type>(invert ? -1 : 0));
359
 
            }
360
 
            else
361
 
            {
362
 
               // For largers values of the non-centrality
363
 
               // parameter Ding's method will consume an
364
 
               // extra-ordinary number of terms, and worse
365
 
               // may return zero when the result is in fact
366
 
               // finite, use Krishnamoorthy's method instead:
367
 
               result = detail::non_central_chi_square_p(
368
 
                  static_cast<value_type>(x), 
369
 
                  static_cast<value_type>(k), 
370
 
                  static_cast<value_type>(l), 
371
 
                  forwarding_policy(),
372
 
                  static_cast<value_type>(invert ? -1 : 0));
373
 
            }
374
 
            if(invert)
375
 
               result = -result;
376
 
            return policies::checked_narrowing_cast<RealType, forwarding_policy>(
377
 
               result, 
378
 
               "boost::math::non_central_chi_squared_cdf<%1%>(%1%, %1%, %1%)");
379
 
         }
380
 
 
381
 
         template <class T, class Policy>
382
 
         struct nccs_quantile_functor
383
 
         {
384
 
            nccs_quantile_functor(const non_central_chi_squared_distribution<T,Policy>& d, T t, bool c)
385
 
               : dist(d), target(t), comp(c) {}
386
 
 
387
 
            T operator()(const T& x)
388
 
            {
389
 
               return comp ?
390
 
                  target - cdf(complement(dist, x))
391
 
                  : cdf(dist, x) - target;
392
 
            }
393
 
 
394
 
         private:
395
 
            non_central_chi_squared_distribution<T,Policy> dist;
396
 
            T target;
397
 
            bool comp;
398
 
         };
399
 
 
400
 
         template <class RealType, class Policy>
401
 
         RealType nccs_quantile(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& p, bool comp)
402
 
         {
403
 
            static const char* function = "quantile(non_central_chi_squared_distribution<%1%>, %1%)";
404
 
            typedef typename policies::evaluation<RealType, Policy>::type value_type;
405
 
            typedef typename policies::normalise<
406
 
               Policy, 
407
 
               policies::promote_float<false>, 
408
 
               policies::promote_double<false>, 
409
 
               policies::discrete_quantile<>,
410
 
               policies::assert_undefined<> >::type forwarding_policy;
411
 
 
412
 
            value_type k = dist.degrees_of_freedom();
413
 
            value_type l = dist.non_centrality();
414
 
            value_type r;
415
 
            if(!detail::check_df(
416
 
               function,
417
 
               k, &r, Policy())
418
 
               ||
419
 
            !detail::check_non_centrality(
420
 
               function,
421
 
               l,
422
 
               &r,
423
 
               Policy())
424
 
               ||
425
 
            !detail::check_probability(
426
 
               function,
427
 
               static_cast<value_type>(p),
428
 
               &r,
429
 
               Policy()))
430
 
                  return (RealType)r;
431
 
 
432
 
            value_type b = (l * l) / (k + 3 * l);
433
 
            value_type c = (k + 3 * l) / (k + 2 * l);
434
 
            value_type ff = (k + 2 * l) / (c * c);
435
 
            value_type guess;
436
 
            if(comp)
437
 
               guess = b + c * quantile(complement(chi_squared_distribution<value_type, forwarding_policy>(ff), p));
438
 
            else
439
 
               guess = b + c * quantile(chi_squared_distribution<value_type, forwarding_policy>(ff), p);
440
 
 
441
 
            if(guess < 0)
442
 
               guess = tools::min_value<value_type>();
443
 
 
444
 
            value_type result = detail::generic_quantile(
445
 
               non_central_chi_squared_distribution<value_type, forwarding_policy>(k, l), 
446
 
               p, 
447
 
               guess, 
448
 
               comp, 
449
 
               function);
450
 
            return policies::checked_narrowing_cast<RealType, forwarding_policy>(
451
 
               result, 
452
 
               function);
453
 
         }
454
 
 
455
 
         template <class RealType, class Policy>
456
 
         RealType nccs_pdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
457
 
         {
458
 
            BOOST_MATH_STD_USING
459
 
            static const char* function = "pdf(non_central_chi_squared_distribution<%1%>, %1%)";
460
 
            typedef typename policies::evaluation<RealType, Policy>::type value_type;
461
 
            typedef typename policies::normalise<
462
 
               Policy, 
463
 
               policies::promote_float<false>, 
464
 
               policies::promote_double<false>, 
465
 
               policies::discrete_quantile<>,
466
 
               policies::assert_undefined<> >::type forwarding_policy;
467
 
 
468
 
            value_type k = dist.degrees_of_freedom();
469
 
            value_type l = dist.non_centrality();
470
 
            value_type r;
471
 
            if(!detail::check_df(
472
 
               function,
473
 
               k, &r, Policy())
474
 
               ||
475
 
            !detail::check_non_centrality(
476
 
               function,
477
 
               l,
478
 
               &r,
479
 
               Policy())
480
 
               ||
481
 
            !detail::check_positive_x(
482
 
               function,
483
 
               (value_type)x,
484
 
               &r,
485
 
               Policy()))
486
 
                  return (RealType)r;
487
 
 
488
 
         BOOST_MATH_STD_USING
489
 
         if(l == 0)
490
 
            return pdf(boost::math::chi_squared_distribution<RealType, forwarding_policy>(dist.degrees_of_freedom()), x);
491
 
 
492
 
         // Special case:
493
 
         if(x == 0)
494
 
            return 0;
495
 
         if(l > 50)
496
 
         {
497
 
            r = non_central_chi_square_pdf(static_cast<value_type>(x), k, l, forwarding_policy());
498
 
         }
499
 
         else
500
 
         {
501
 
            r = log(x / l) * (k / 4 - 0.5f) - (x + l) / 2;
502
 
            if(fabs(r) >= tools::log_max_value<RealType>() / 4)
503
 
            {
504
 
               r = non_central_chi_square_pdf(static_cast<value_type>(x), k, l, forwarding_policy());
505
 
            }
506
 
            else
507
 
            {
508
 
               r = exp(r);
509
 
               r = 0.5f * r
510
 
                  * boost::math::cyl_bessel_i(k/2 - 1, sqrt(l * x), forwarding_policy());
511
 
            }
512
 
         }
513
 
         return policies::checked_narrowing_cast<RealType, forwarding_policy>(
514
 
               r, 
515
 
               function);
516
 
         }
517
 
 
518
 
         template <class RealType, class Policy>
519
 
         struct degrees_of_freedom_finder
520
 
         {
521
 
            degrees_of_freedom_finder(
522
 
               RealType lam_, RealType x_, RealType p_, bool c)
523
 
               : lam(lam_), x(x_), p(p_), comp(c) {}
524
 
 
525
 
            RealType operator()(const RealType& v)
526
 
            {
527
 
               non_central_chi_squared_distribution<RealType, Policy> d(v, lam);
528
 
               return comp ?
529
 
                  p - cdf(complement(d, x))
530
 
                  : cdf(d, x) - p;
531
 
            }
532
 
         private:
533
 
            RealType lam;
534
 
            RealType x;
535
 
            RealType p;
536
 
            bool comp;
537
 
         };
538
 
 
539
 
         template <class RealType, class Policy>
540
 
         inline RealType find_degrees_of_freedom(
541
 
            RealType lam, RealType x, RealType p, RealType q, const Policy& pol)
542
 
         {
543
 
            const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom";
544
 
            if((p == 0) || (q == 0))
545
 
            {
546
 
               //
547
 
               // Can't a thing if one of p and q is zero:
548
 
               //
549
 
               return policies::raise_evaluation_error<RealType>(function, 
550
 
                  "Can't find degrees of freedom when the probability is 0 or 1, only possible answer is %1%", 
551
 
                  RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy());
552
 
            }
553
 
            degrees_of_freedom_finder<RealType, Policy> f(lam, x, p < q ? p : q, p < q ? false : true);
554
 
            tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
555
 
            boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
556
 
            //
557
 
            // Pick an initial guess that we know will give us a probability
558
 
            // right around 0.5.
559
 
            //
560
 
            RealType guess = x - lam;
561
 
            if(guess < 1)
562
 
               guess = 1;
563
 
            std::pair<RealType, RealType> ir = tools::bracket_and_solve_root(
564
 
               f, guess, RealType(2), false, tol, max_iter, pol);
565
 
            RealType result = ir.first + (ir.second - ir.first) / 2;
566
 
            if(max_iter >= policies::get_max_root_iterations<Policy>())
567
 
            {
568
 
               policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
569
 
                  " or there is no answer to problem.  Current best guess is %1%", result, Policy());
570
 
            }
571
 
            return result;
572
 
         }
573
 
 
574
 
         template <class RealType, class Policy>
575
 
         struct non_centrality_finder
576
 
         {
577
 
            non_centrality_finder(
578
 
               RealType v_, RealType x_, RealType p_, bool c)
579
 
               : v(v_), x(x_), p(p_), comp(c) {}
580
 
 
581
 
            RealType operator()(const RealType& lam)
582
 
            {
583
 
               non_central_chi_squared_distribution<RealType, Policy> d(v, lam);
584
 
               return comp ?
585
 
                  p - cdf(complement(d, x))
586
 
                  : cdf(d, x) - p;
587
 
            }
588
 
         private:
589
 
            RealType v;
590
 
            RealType x;
591
 
            RealType p;
592
 
            bool comp;
593
 
         };
594
 
 
595
 
         template <class RealType, class Policy>
596
 
         inline RealType find_non_centrality(
597
 
            RealType v, RealType x, RealType p, RealType q, const Policy& pol)
598
 
         {
599
 
            const char* function = "non_central_chi_squared<%1%>::find_non_centrality";
600
 
            if((p == 0) || (q == 0))
601
 
            {
602
 
               //
603
 
               // Can't do a thing if one of p and q is zero:
604
 
               //
605
 
               return policies::raise_evaluation_error<RealType>(function, 
606
 
                  "Can't find non centrality parameter when the probability is 0 or 1, only possible answer is %1%", 
607
 
                  RealType(std::numeric_limits<RealType>::quiet_NaN()), Policy());
608
 
            }
609
 
            non_centrality_finder<RealType, Policy> f(v, x, p < q ? p : q, p < q ? false : true);
610
 
            tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
611
 
            boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
612
 
            //
613
 
            // Pick an initial guess that we know will give us a probability
614
 
            // right around 0.5.
615
 
            //
616
 
            RealType guess = x - v;
617
 
            if(guess < 1)
618
 
               guess = 1;
619
 
            std::pair<RealType, RealType> ir = tools::bracket_and_solve_root(
620
 
               f, guess, RealType(2), false, tol, max_iter, pol);
621
 
            RealType result = ir.first + (ir.second - ir.first) / 2;
622
 
            if(max_iter >= policies::get_max_root_iterations<Policy>())
623
 
            {
624
 
               policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
625
 
                  " or there is no answer to problem.  Current best guess is %1%", result, Policy());
626
 
            }
627
 
            return result;
628
 
         }
629
 
 
630
 
      }
631
 
 
632
 
      template <class RealType = double, class Policy = policies::policy<> >
633
 
      class non_central_chi_squared_distribution
634
 
      {
635
 
      public:
636
 
         typedef RealType value_type;
637
 
         typedef Policy policy_type;
638
 
 
639
 
         non_central_chi_squared_distribution(RealType df_, RealType lambda) : df(df_), ncp(lambda)
640
 
         { 
641
 
            const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::non_central_chi_squared_distribution(%1%,%1%)";
642
 
            RealType r;
643
 
            detail::check_df(
644
 
               function,
645
 
               df, &r, Policy());
646
 
            detail::check_non_centrality(
647
 
               function,
648
 
               ncp,
649
 
               &r,
650
 
               Policy());
651
 
         } // non_central_chi_squared_distribution constructor.
652
 
 
653
 
         RealType degrees_of_freedom() const
654
 
         { // Private data getter function.
655
 
            return df;
656
 
         }
657
 
         RealType non_centrality() const
658
 
         { // Private data getter function.
659
 
            return ncp;
660
 
         }
661
 
         static RealType find_degrees_of_freedom(RealType lam, RealType x, RealType p)
662
 
         {
663
 
            const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom";
664
 
            typedef typename policies::evaluation<RealType, Policy>::type value_type;
665
 
            typedef typename policies::normalise<
666
 
               Policy, 
667
 
               policies::promote_float<false>, 
668
 
               policies::promote_double<false>, 
669
 
               policies::discrete_quantile<>,
670
 
               policies::assert_undefined<> >::type forwarding_policy;
671
 
            value_type result = detail::find_degrees_of_freedom(
672
 
               static_cast<value_type>(lam), 
673
 
               static_cast<value_type>(x), 
674
 
               static_cast<value_type>(p), 
675
 
               static_cast<value_type>(1-p), 
676
 
               forwarding_policy());
677
 
            return policies::checked_narrowing_cast<RealType, forwarding_policy>(
678
 
               result, 
679
 
               function);
680
 
         }
681
 
         template <class A, class B, class C>
682
 
         static RealType find_degrees_of_freedom(const complemented3_type<A,B,C>& c)
683
 
         {
684
 
            const char* function = "non_central_chi_squared<%1%>::find_degrees_of_freedom";
685
 
            typedef typename policies::evaluation<RealType, Policy>::type value_type;
686
 
            typedef typename policies::normalise<
687
 
               Policy, 
688
 
               policies::promote_float<false>, 
689
 
               policies::promote_double<false>, 
690
 
               policies::discrete_quantile<>,
691
 
               policies::assert_undefined<> >::type forwarding_policy;
692
 
            value_type result = detail::find_degrees_of_freedom(
693
 
               static_cast<value_type>(c.dist), 
694
 
               static_cast<value_type>(c.param1), 
695
 
               static_cast<value_type>(1-c.param2), 
696
 
               static_cast<value_type>(c.param2), 
697
 
               forwarding_policy());
698
 
            return policies::checked_narrowing_cast<RealType, forwarding_policy>(
699
 
               result, 
700
 
               function);
701
 
         }
702
 
         static RealType find_non_centrality(RealType v, RealType x, RealType p)
703
 
         {
704
 
            const char* function = "non_central_chi_squared<%1%>::find_non_centrality";
705
 
            typedef typename policies::evaluation<RealType, Policy>::type value_type;
706
 
            typedef typename policies::normalise<
707
 
               Policy, 
708
 
               policies::promote_float<false>, 
709
 
               policies::promote_double<false>, 
710
 
               policies::discrete_quantile<>,
711
 
               policies::assert_undefined<> >::type forwarding_policy;
712
 
            value_type result = detail::find_non_centrality(
713
 
               static_cast<value_type>(v), 
714
 
               static_cast<value_type>(x), 
715
 
               static_cast<value_type>(p), 
716
 
               static_cast<value_type>(1-p), 
717
 
               forwarding_policy());
718
 
            return policies::checked_narrowing_cast<RealType, forwarding_policy>(
719
 
               result, 
720
 
               function);
721
 
         }
722
 
         template <class A, class B, class C>
723
 
         static RealType find_non_centrality(const complemented3_type<A,B,C>& c)
724
 
         {
725
 
            const char* function = "non_central_chi_squared<%1%>::find_non_centrality";
726
 
            typedef typename policies::evaluation<RealType, Policy>::type value_type;
727
 
            typedef typename policies::normalise<
728
 
               Policy, 
729
 
               policies::promote_float<false>, 
730
 
               policies::promote_double<false>, 
731
 
               policies::discrete_quantile<>,
732
 
               policies::assert_undefined<> >::type forwarding_policy;
733
 
            value_type result = detail::find_non_centrality(
734
 
               static_cast<value_type>(c.dist), 
735
 
               static_cast<value_type>(c.param1), 
736
 
               static_cast<value_type>(1-c.param2), 
737
 
               static_cast<value_type>(c.param2), 
738
 
               forwarding_policy());
739
 
            return policies::checked_narrowing_cast<RealType, forwarding_policy>(
740
 
               result, 
741
 
               function);
742
 
         }
743
 
      private:
744
 
         // Data member, initialized by constructor.
745
 
         RealType df; // degrees of freedom.
746
 
         RealType ncp; // non-centrality parameter
747
 
      }; // template <class RealType, class Policy> class non_central_chi_squared_distribution
748
 
 
749
 
      typedef non_central_chi_squared_distribution<double> non_central_chi_squared; // Reserved name of type double.
750
 
 
751
 
      // Non-member functions to give properties of the distribution.
752
 
 
753
 
      template <class RealType, class Policy>
754
 
      inline const std::pair<RealType, RealType> range(const non_central_chi_squared_distribution<RealType, Policy>& /* dist */)
755
 
      { // Range of permissible values for random variable k.
756
 
         using boost::math::tools::max_value;
757
 
         return std::pair<RealType, RealType>(0, max_value<RealType>()); // Max integer?
758
 
      }
759
 
 
760
 
      template <class RealType, class Policy>
761
 
      inline const std::pair<RealType, RealType> support(const non_central_chi_squared_distribution<RealType, Policy>& /* dist */)
762
 
      { // Range of supported values for random variable k.
763
 
         // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
764
 
         using boost::math::tools::max_value;
765
 
         return std::pair<RealType, RealType>(0,  max_value<RealType>());
766
 
      }
767
 
 
768
 
      template <class RealType, class Policy>
769
 
      inline RealType mean(const non_central_chi_squared_distribution<RealType, Policy>& dist)
770
 
      { // Mean of poisson distribution = lambda.
771
 
         const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::mean()";
772
 
         RealType k = dist.degrees_of_freedom();
773
 
         RealType l = dist.non_centrality();
774
 
         RealType r;
775
 
         if(!detail::check_df(
776
 
            function,
777
 
            k, &r, Policy())
778
 
            ||
779
 
         !detail::check_non_centrality(
780
 
            function,
781
 
            l,
782
 
            &r,
783
 
            Policy()))
784
 
               return r;
785
 
         return k + l;
786
 
      } // mean
787
 
 
788
 
      template <class RealType, class Policy>
789
 
      inline RealType mode(const non_central_chi_squared_distribution<RealType, Policy>& dist)
790
 
      { // mode.
791
 
         static const char* function = "mode(non_central_chi_squared_distribution<%1%> const&)";
792
 
 
793
 
         RealType k = dist.degrees_of_freedom();
794
 
         RealType l = dist.non_centrality();
795
 
         RealType r;
796
 
         if(!detail::check_df(
797
 
            function,
798
 
            k, &r, Policy())
799
 
            ||
800
 
         !detail::check_non_centrality(
801
 
            function,
802
 
            l,
803
 
            &r,
804
 
            Policy()))
805
 
               return (RealType)r;
806
 
         return detail::generic_find_mode(dist, 1 + k, function);
807
 
      }
808
 
 
809
 
      template <class RealType, class Policy>
810
 
      inline RealType variance(const non_central_chi_squared_distribution<RealType, Policy>& dist)
811
 
      { // variance.
812
 
         const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::variance()";
813
 
         RealType k = dist.degrees_of_freedom();
814
 
         RealType l = dist.non_centrality();
815
 
         RealType r;
816
 
         if(!detail::check_df(
817
 
            function,
818
 
            k, &r, Policy())
819
 
            ||
820
 
         !detail::check_non_centrality(
821
 
            function,
822
 
            l,
823
 
            &r,
824
 
            Policy()))
825
 
               return r;
826
 
         return 2 * (2 * l + k);
827
 
      }
828
 
 
829
 
      // RealType standard_deviation(const non_central_chi_squared_distribution<RealType, Policy>& dist)
830
 
      // standard_deviation provided by derived accessors.
831
 
 
832
 
      template <class RealType, class Policy>
833
 
      inline RealType skewness(const non_central_chi_squared_distribution<RealType, Policy>& dist)
834
 
      { // skewness = sqrt(l).
835
 
         const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::skewness()";
836
 
         RealType k = dist.degrees_of_freedom();
837
 
         RealType l = dist.non_centrality();
838
 
         RealType r;
839
 
         if(!detail::check_df(
840
 
            function,
841
 
            k, &r, Policy())
842
 
            ||
843
 
         !detail::check_non_centrality(
844
 
            function,
845
 
            l,
846
 
            &r,
847
 
            Policy()))
848
 
               return r;
849
 
         BOOST_MATH_STD_USING
850
 
            return pow(2 / (k + 2 * l), RealType(3)/2) * (k + 3 * l);
851
 
      }
852
 
 
853
 
      template <class RealType, class Policy>
854
 
      inline RealType kurtosis_excess(const non_central_chi_squared_distribution<RealType, Policy>& dist)
855
 
      { 
856
 
         const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::kurtosis_excess()";
857
 
         RealType k = dist.degrees_of_freedom();
858
 
         RealType l = dist.non_centrality();
859
 
         RealType r;
860
 
         if(!detail::check_df(
861
 
            function,
862
 
            k, &r, Policy())
863
 
            ||
864
 
         !detail::check_non_centrality(
865
 
            function,
866
 
            l,
867
 
            &r,
868
 
            Policy()))
869
 
               return r;
870
 
         return 12 * (k + 4 * l) / ((k + 2 * l) * (k + 2 * l));
871
 
      } // kurtosis_excess
872
 
 
873
 
      template <class RealType, class Policy>
874
 
      inline RealType kurtosis(const non_central_chi_squared_distribution<RealType, Policy>& dist)
875
 
      {
876
 
         return kurtosis_excess(dist) + 3;
877
 
      }
878
 
 
879
 
      template <class RealType, class Policy>
880
 
      inline RealType pdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
881
 
      { // Probability Density/Mass Function.
882
 
         return detail::nccs_pdf(dist, x);
883
 
      } // pdf
884
 
 
885
 
      template <class RealType, class Policy>
886
 
      RealType cdf(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
887
 
      { 
888
 
         const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::cdf(%1%)";
889
 
         RealType k = dist.degrees_of_freedom();
890
 
         RealType l = dist.non_centrality();
891
 
         RealType r;
892
 
         if(!detail::check_df(
893
 
            function,
894
 
            k, &r, Policy())
895
 
            ||
896
 
         !detail::check_non_centrality(
897
 
            function,
898
 
            l,
899
 
            &r,
900
 
            Policy())
901
 
            ||
902
 
         !detail::check_positive_x(
903
 
            function,
904
 
            x,
905
 
            &r,
906
 
            Policy()))
907
 
               return r;
908
 
 
909
 
         return detail::non_central_chi_squared_cdf(x, k, l, false, Policy());
910
 
      } // cdf
911
 
 
912
 
      template <class RealType, class Policy>
913
 
      RealType cdf(const complemented2_type<non_central_chi_squared_distribution<RealType, Policy>, RealType>& c)
914
 
      { // Complemented Cumulative Distribution Function
915
 
         const char* function = "boost::math::non_central_chi_squared_distribution<%1%>::cdf(%1%)";
916
 
         non_central_chi_squared_distribution<RealType, Policy> const& dist = c.dist;
917
 
         RealType x = c.param;
918
 
         RealType k = dist.degrees_of_freedom();
919
 
         RealType l = dist.non_centrality();
920
 
         RealType r;
921
 
         if(!detail::check_df(
922
 
            function,
923
 
            k, &r, Policy())
924
 
            ||
925
 
         !detail::check_non_centrality(
926
 
            function,
927
 
            l,
928
 
            &r,
929
 
            Policy())
930
 
            ||
931
 
         !detail::check_positive_x(
932
 
            function,
933
 
            x,
934
 
            &r,
935
 
            Policy()))
936
 
               return r;
937
 
 
938
 
         return detail::non_central_chi_squared_cdf(x, k, l, true, Policy());
939
 
      } // ccdf
940
 
 
941
 
      template <class RealType, class Policy>
942
 
      inline RealType quantile(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& p)
943
 
      { // Quantile (or Percent Point) function.
944
 
         return detail::nccs_quantile(dist, p, false);
945
 
      } // quantile
946
 
 
947
 
      template <class RealType, class Policy>
948
 
      inline RealType quantile(const complemented2_type<non_central_chi_squared_distribution<RealType, Policy>, RealType>& c)
949
 
      { // Quantile (or Percent Point) function.
950
 
         return detail::nccs_quantile(c.dist, c.param, true);
951
 
      } // quantile complement.
952
 
 
953
 
   } // namespace math
954
 
} // namespace boost
955
 
 
956
 
// This include must be at the end, *after* the accessors
957
 
// for this distribution have been defined, in order to
958
 
// keep compilers that support two-phase lookup happy.
959
 
#include <boost/math/distributions/detail/derived_accessors.hpp>
960
 
 
961
 
#endif // BOOST_MATH_SPECIAL_NON_CENTRAL_CHI_SQUARE_HPP
962
 
 
963
 
 
964