chisquaredistribution.cpp

有很多的函数库

/* -*- mode: c++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */

/*
 Copyright (C) 2002, 2003 Sadruddin Rejeb
 Copyright (C) 2007 Klaus Spanderen

 This file is part of QuantLib, a free-software/open-source library
 for financial quantitative analysts and developers - http://quantlib.org/

 QuantLib is free software: you can redistribute it and/or modify it
 under the terms of the QuantLib license.  You should have received a
 copy of the license along with this program; if not, please email
 <[e-mail]>. The license is also available online at
 <http://quantlib.org/license.shtml>.

 This program is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
 FOR A PARTICULAR PURPOSE.  See the license for more details.
*/

#include <ql/math/solvers1d/brent.hpp>
#include <ql/math/functional.hpp>
#include <ql/math/distributions/chisquaredistribution.hpp>
#include <ql/math/distributions/gammadistribution.hpp>
#include <ql/math/distributions/normaldistribution.hpp>

namespace QuantLib {

    Real ChiSquareDistribution::operator()(Real x) const {
        return GammaDistribution(0.5*df_)(0.5*x);
    }

    Real NonCentralChiSquareDistribution::operator()(Real x) const {
        if (x <= 0.0)
            return 0.0;

        const Real errmax = 1e-12;
        const Size itrmax = 10000;
        Real lam = 0.5*ncp_;

        Real u = std::exp(-lam);
        Real v = u;
        Real x2 = 0.5*x;
        Real f2 = 0.5*df_;
        Real f_x_2n = df_ - x;

        Real t = 0.0;
        if (f2*QL_EPSILON > 0.125 &&
            std::fabs(x2-f2) < std::sqrt(QL_EPSILON)*f2) {
            t = std::exp((1 - t) *
                         (2 - t/(f2+1)))/std::sqrt(2.0*M_PI*(f2 + 1.0));
        }
        else {
            t = std::exp(f2*std::log(x2) - x2 -
                         GammaFunction().logValue(f2 + 1));
        }

        Real ans = v*t;

        bool flag = false;
        Size n = 1;
        Real f_2n = df_ + 2.0;
        f_x_2n += 2.0;

        Real bound;
        for (;;) {
            if (f_x_2n > 0) {
                flag = true;
                goto L10;
            }
            for (;;) {
                u *= lam / n;
                v += u;
                t *= x / f_2n;
                ans += v*t;
                n++;
                f_2n += 2.0;
                f_x_2n += 2.0;
                if (!flag && n <= itrmax)
                    break;
            L10:
                bound = t * x / f_x_2n;
                if (bound <= errmax || n > itrmax)
                    goto L_End;
            }
        }
    L_End:
        if (bound > errmax) QL_FAIL("didn't converge");
        return (ans);

    }

    InverseNonCentralChiSquareDistribution::
      InverseNonCentralChiSquareDistribution(Real df, Real ncp, 
                                             Size maxEvaluations, 
                                             Real accuracy)
    : nonCentralDist_(df, ncp),
      guess_(df+ncp),
      maxEvaluations_(maxEvaluations),
      accuracy_(accuracy) {
    }

    Real InverseNonCentralChiSquareDistribution::operator()(Real x) const {

        // first find the right side of the interval
        Real upper = guess_;
        Size evaluations = maxEvaluations_;
        while (nonCentralDist_(upper) < x && evaluations > 0) {
            upper*=2.0;
            --evaluations;
        }

        // use a brent solver for the rest
        Brent solver;
        solver.setMaxEvaluations(evaluations);
        return solver.solve(compose(std::bind2nd(std::minus<Real>(),x), 
                                    nonCentralDist_),
                            accuracy_, 0.75*upper, 
                            (evaluations == maxEvaluations_)? 0.0: 0.5*upper,
                            upper);
    }
}