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 27 namespace QuantLib { 28 operator ()(Real x) const29 Real CumulativeChiSquareDistribution::operator()(Real x) const { 30 return CumulativeGammaDistribution(0.5*df_)(0.5*x); 31 } 32 operator ()(Real x) const33 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(-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(x2-f2) < std::sqrt(QL_EPSILON)*f2) { 50 t = std::exp((1 - t) * 51 (2 - t/(f2+1)))/std::sqrt(2.0*M_PI*(f2 + 1.0)); 52 } 53 else { 54 t = std::exp(f2*std::log(x2) - x2 - 55 GammaFunction().logValue(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 operator ()(Real x) const93 Real NonCentralCumulativeChiSquareSankaranApprox::operator()(Real x) const { 94 95 const Real h = 1-2*(df_+ncp_)*(df_+3*ncp_)/(3*square<Real>()(df_+2*ncp_)); 96 const Real p = (df_+2*ncp_)/square<Real>()(df_+ncp_); 97 const Real m = (h-1)*(1-3*h); 98 99 const Real u= (std::pow(x/(df_+ncp_), h) - (1 + h*p*(h-1-0.5*(2-h)*m*p)))/ 100 (h*std::sqrt(2*p)*(1+0.5*m*p)); 101 102 return CumulativeNormalDistribution()(u); 103 } 104 105 InverseNonCentralCumulativeChiSquareDistribution:: InverseNonCentralCumulativeChiSquareDistribution(Real df,Real ncp,Size maxEvaluations,Real accuracy)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 operator ()(Real x) const115 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(compose(subtract<Real>(x), 129 nonCentralDist_), 130 accuracy_, 0.75*upper, 131 (evaluations == maxEvaluations_)? 0.0: 0.5*upper, 132 upper); 133 } 134 } 135