1/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2
3/*
4 Copyright (C) 2002, 2003 Sadruddin Rejeb
5 Copyright (C) 2007 Klaus Spanderen
6
7 This file is part of QuantLib, a free-software/open-source library
8 for financial quantitative analysts and developers - http://quantlib.org/
9
10 QuantLib is free software: you can redistribute it and/or modify it
11 under the terms of the QuantLib license. You should have received a
12 copy of the license along with this program; if not, please email
13 <quantlib-dev@lists.sf.net>. The license is also available online at
14 <http://quantlib.org/license.shtml>.
15
16 This program is distributed in the hope that it will be useful, but WITHOUT
17 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 FOR A PARTICULAR PURPOSE. See the license for more details.
19*/
20
21#include <ql/math/solvers1d/brent.hpp>
22#include <ql/math/functional.hpp>
23#include <ql/math/distributions/chisquaredistribution.hpp>
24#include <ql/math/distributions/gammadistribution.hpp>
25#include <ql/math/distributions/normaldistribution.hpp>
26
27namespace QuantLib {
28
29 Real CumulativeChiSquareDistribution::operator()(Real x) const {
30 return CumulativeGammaDistribution(0.5*df_)(0.5*x);
31 }
32
33 Real NonCentralCumulativeChiSquareDistribution::operator()(Real x) const {
34 if (x <= 0.0)
35 return 0.0;
36
37 const Real errmax = 1e-12;
38 const Size itrmax = 10000;
39 Real lam = 0.5*ncp_;
40
41 Real u = std::exp(x: -lam);
42 Real v = u;
43 Real x2 = 0.5*x;
44 Real f2 = 0.5*df_;
45 Real f_x_2n = df_ - x;
46
47 Real t = 0.0;
48 if (f2*QL_EPSILON > 0.125 &&
49 std::fabs(x: x2-f2) < std::sqrt(QL_EPSILON)*f2) {
50 t = std::exp(x: (1 - t) *
51 (2 - t/(f2+1)))/std::sqrt(x: 2.0*M_PI*(f2 + 1.0));
52 }
53 else {
54 t = std::exp(x: f2*std::log(x: x2) - x2 -
55 GammaFunction().logValue(x: f2 + 1));
56 }
57
58 Real ans = v*t;
59
60 bool flag = false;
61 Size n = 1;
62 Real f_2n = df_ + 2.0;
63 f_x_2n += 2.0;
64
65 Real bound;
66 for (;;) {
67 if (f_x_2n > 0) {
68 flag = true;
69 goto L10;
70 }
71 for (;;) {
72 u *= lam / n;
73 v += u;
74 t *= x / f_2n;
75 ans += v*t;
76 n++;
77 f_2n += 2.0;
78 f_x_2n += 2.0;
79 if (!flag && n <= itrmax)
80 break;
81 L10:
82 bound = t * x / f_x_2n;
83 if (bound <= errmax || n > itrmax)
84 goto L_End;
85 }
86 }
87 L_End:
88 if (bound > errmax) QL_FAIL("didn't converge");
89 return (ans);
90
91 }
92
93 Real NonCentralCumulativeChiSquareSankaranApprox::operator()(Real x) const {
94
95 const Real h = 1-2*(df_+ncp_)*(df_+3*ncp_)/(3*squared(x: df_+2*ncp_));
96 const Real p = (df_+2*ncp_)/squared(x: df_+ncp_);
97 const Real m = (h-1)*(1-3*h);
98
99 const Real u= (std::pow(x: x/(df_+ncp_), y: h) - (1 + h*p*(h-1-0.5*(2-h)*m*p)))/
100 (h*std::sqrt(x: 2*p)*(1+0.5*m*p));
101
102 return CumulativeNormalDistribution()(u);
103 }
104
105 InverseNonCentralCumulativeChiSquareDistribution::
106 InverseNonCentralCumulativeChiSquareDistribution(Real df, Real ncp,
107 Size maxEvaluations,
108 Real accuracy)
109 : nonCentralDist_(df, ncp),
110 guess_(df+ncp),
111 maxEvaluations_(maxEvaluations),
112 accuracy_(accuracy) {
113 }
114
115 Real InverseNonCentralCumulativeChiSquareDistribution::operator()(Real x) const {
116
117 // first find the right side of the interval
118 Real upper = guess_;
119 Size evaluations = maxEvaluations_;
120 while (nonCentralDist_(upper) < x && evaluations > 0) {
121 upper*=2.0;
122 --evaluations;
123 }
124
125 // use a Brent solver for the rest
126 Brent solver;
127 solver.setMaxEvaluations(evaluations);
128 return solver.solve(f: [&](Real y) { return nonCentralDist_(y) - x; },
129 accuracy: accuracy_, guess: 0.75*upper,
130 xMin: (evaluations == maxEvaluations_)? 0.0: Real(0.5*upper),
131 xMax: upper);
132 }
133}
134

source code of quantlib/ql/math/distributions/chisquaredistribution.cpp