| /* ----------------------------------------------------------------------- *//** |
| * |
| * @file student.hpp |
| * |
| * @brief Evaluate the Student's t-distribution function. |
| * @author Florian Schoppmann |
| * @date November 2010 |
| * |
| *//* -------------------------------------------------------------------- *//** |
| * |
| * @file student.hpp |
| * |
| * Emprirical results indicate that the numerical quality of the series |
| * expansion from [1] (see notes below) is vastly superior to using continued |
| * fractions for computing the cdf via the incomplete beta function. |
| * |
| * @literature |
| * |
| * [1] Abramowitz and Stegun, Handbook of Mathematical Functions with Formulas, |
| * Graphs, and Mathematical Tables, 1972 |
| * page 948: http://people.math.sfu.ca/~cbm/aands/page_948.htm |
| * |
| * Further reading (for computing the Student-T cdf via the incomplete beta |
| * function): |
| * |
| * [2] NIST Digital Library of Mathematical Functions, Ch. 8, |
| * Incomplete Gamma and Related Functions, |
| * http://dlmf.nist.gov/8.17 |
| * |
| * [3] Lentz, Generating Bessel functions in Mie scattering calculations using |
| * continued fractions, Applied Optics, Vol. 15, No. 3, 1976 |
| * |
| * [4] Thompson and Barnett, Coulomb and Bessel Functions of Complex Arguments |
| * and Order, Journal of Computational Physics, Vol. 64, 1986 |
| * |
| * [5] Cuyt et al., Handbook of Continued Fractions for Special Functions, |
| * Springer, 2008 |
| * |
| * [6] Gil et al., Numerical Methods for Special Functions, SIAM, 2008 |
| * |
| * [7] Press et al., Numerical Recipes in C++, 3rd edition, |
| * Cambridge Univ. Press, 2007 |
| * |
| * [8] DiDonato, Morris, Jr., Algorithm 708: Significant Digit Computation of |
| * the Incomplete Beta Function Ratios, ACM Transactions on Mathematical |
| * Software, Vol. 18, No. 3, 1992 |
| * |
| * Approximating the Student-T distribution function with the normal |
| * distribution: |
| * |
| * [9] Gleason, A note on a proposed student t approximation, Computational |
| * Statistics & Data Analysis, Vol. 34, No. 1, 2000 |
| * |
| * [10] Gaver and Kafadar, A Retrievable Recipe for Inverse t, The American |
| * Statistician, Vol. 38, No. 4, 1984 |
| */ |
| |
| /** |
| * @brief Student-t cumulative distribution function |
| */ |
| DECLARE_UDF(prob, students_t_cdf) |
| DECLARE_UDF(prob, students_t_pdf) |
| DECLARE_UDF(prob, students_t_quantile) |
| |
| |
| #ifndef MADLIB_MODULES_PROB_STUDENT_T_HPP |
| #define MADLIB_MODULES_PROB_STUDENT_T_HPP |
| |
| #include <boost/math/distributions/detail/common_error_handling.hpp> |
| #include <boost/math/distributions/normal.hpp> |
| #include <boost/math/distributions/students_t.hpp> |
| |
| namespace madlib { |
| |
| namespace modules { |
| |
| namespace prob { |
| |
| typedef boost::math::students_t_distribution<double, boost_mathkit_policy> |
| students_t; |
| |
| namespace { |
| |
| /** |
| * @brief Compute one-sided Student's t cumulative distribution function |
| * |
| * We use the series expansions 26.7.3 and 26.7.4 from [1] and |
| * substitute sin(theta) = t/sqrt(n * z), where z = 1 + t^2/nu. |
| * |
| * This gives: |
| * @verbatim |
| * t |
| * A(t|1) = 2 arctan( -------- ) , |
| * sqrt(nu) |
| * |
| * (nu-3)/2 |
| * 2 [ t t -- 2 * 4 * ... * (2i) ] |
| * A(t|nu) = - * [ arctan( -------- ) + ------------ * \ ---------------------- ] |
| * π [ sqrt(nu) sqrt(nu) * z /_ 3 * ... * (2i+1) * z^i ] |
| * i=0 |
| * for odd nu > 1, and |
| * |
| * (nu-2)/2 |
| * t -- 1 * 3 * ... * (2i - 1) |
| * A(t|nu) = ------------ * \ ------------------------ for even nu, |
| * sqrt(nu * z) /_ 2 * 4 * ... * (2i) * z^i |
| * i=0 |
| * |
| * where A(t|nu) = Pr[|T| <= t]. |
| * @endverbatim |
| * |
| * @param t |
| * @param nu Degree of freedom \f$ \nu > 0 \f$ |
| * @return \f$ \Pr[|T| < t] \f$ where \f$ t \geq 0 \f$, \f$ T \f$ is a Student's |
| * T-distributed random variable with \f$ \nu \f$ degrees of |
| * freedom. |
| * |
| * Note: The running time of calculating the series is proportional to nu. |
| * We therefore use the normal distribution as an approximation for large nu. |
| * Another idea for handling this case can be found in reference [8]. |
| */ |
| template <class RealType> |
| inline |
| RealType |
| oneSidedStudentsT_CDF(const RealType& t, uint64_t nu) { |
| RealType z, |
| t_by_sqrt_nu; |
| RealType A, /* contains A(t|nu) */ |
| prod = 1., |
| sum = 1.; |
| |
| /* Handle main case (nu \in {1, ..., 200}) in the rest of the function. */ |
| z = 1. + t * t / static_cast<double>(nu); |
| t_by_sqrt_nu = std::fabs(t) / std::sqrt(static_cast<double>(nu)); |
| |
| if (nu == 1) |
| { |
| A = 2. / M_PI * std::atan(t_by_sqrt_nu); |
| } |
| else if (nu & 1) /* odd nu > 1 */ |
| { |
| for (uint64_t j = 2; j + 3 <= nu; j += 2) |
| { |
| prod = prod * static_cast<double>(j) |
| / (static_cast<double>(j + 1) * z); |
| sum = sum + prod; |
| } |
| A = 2 / M_PI * ( std::atan(t_by_sqrt_nu) + t_by_sqrt_nu / z * sum ); |
| } |
| else /* even nu */ |
| { |
| for (uint64_t j = 2; j + 2 <= nu; j += 2) |
| { |
| prod = prod * static_cast<double>(j - 1) |
| / (static_cast<double>(j) * z); |
| sum = sum + prod; |
| } |
| A = t_by_sqrt_nu / std::sqrt(z) * sum; |
| } |
| |
| /* A should obviously be within the interval [0,1] plus minus (hopefully |
| * small) rounding errors. */ |
| if (A > 1.) |
| A = 1.; |
| else if (A < 0.) |
| A = 0.; |
| |
| return A; |
| } |
| |
| /** |
| * @brief Compute parameter for normal CDF for approximating the Student's T CDF |
| * |
| * Gleason suggested a formula for approximating the Student's |
| * t-distribution [9], which goes back to an approximation suggested in [10]. |
| * |
| * Compared to the series expansion, this approximation satisfies |
| * rel_error < 0.0001 || abs_error < 0.00000001 |
| * for all nu >= 200. (Tested on Mac OS X 10.6, gcc-4.2.) |
| * |
| * @param t |
| * @param nu Degree of freedom \f$ \nu > 0 \f$ |
| * @returns A value \f$ z \f$ such that for a Student's t-distributed |
| * random variable \f$ T \f$ with \f$ nu \f$ degrees of freedom and a |
| * standard normally distributed random variable \f$ Z \f$, it holds that |
| * \f$ \Pr[T \leq t] \approx \Pr[Z \leq z] \f$. |
| */ |
| template <class RealType> |
| inline |
| RealType |
| GleasonsNormalApproxForStudentsT(const RealType& t, const RealType& nu) { |
| double g = (nu - 1.5) / ((nu - 1) * (nu - 1)), |
| z = std::sqrt( std::log(1. + t * t / nu) / g ); |
| |
| if (t < 0) |
| z *= -1.; |
| |
| return z; |
| } |
| |
| } // anonymous namespace |
| |
| /** |
| * @brief Compute Student's cumulative distribution function |
| * |
| * For nu >= 1000000, we just use the normal distribution as an approximation. |
| * For 1000000 >= nu >= 200, we use a simple approximation from [9]. |
| * If nu is not within 0.01 of a natural number, we will call the student-t |
| * CDF from boost. Otherwise, our approach should be much more precise than |
| * using the incomplete beta function as boost does (see the references). |
| * |
| * We are much more cautious than usual here (it is folklore that the normal |
| * distribution is a "good" estimate for Student-T if nu >= 30), but we can |
| * afford the extra work as this function is not designed to be called from |
| * inner loops. Performance should still be reasonably good, with at most ~100 |
| * iterations in any case (just one if nu >= 200). |
| * |
| * For nu < 200, we use the series expansions 26.7.3 and 26.7.4 from [1] and |
| * substitute sin(theta) = t/sqrt(n * z), where z = 1 + t^2/nu (using |
| * oneSidedStudentsT_CDF()). |
| * |
| * @param dist A Student's t-distribution object, containing the degree of |
| * freedom \f$ \nu \f$ |
| * @param t |
| * @return \f$ \Pr[T < t] \f$ where \f$ t \geq 0 \f$, \f$ T \f$ is a Student's |
| * T-distributed random variable with \f$ \nu \f$ degrees of |
| * freedom. |
| */ |
| template <class RealType, class Policy> |
| inline |
| RealType |
| cdf(const boost::math::students_t_distribution<RealType, Policy>& dist, |
| const RealType& t) { |
| |
| RealType df = dist.degrees_of_freedom(); |
| |
| // FIXME: Add some justification/do some tests. |
| if (!std::isfinite(df) || std::fabs(df - std::floor(df))/df > 0.01) |
| return boost::math::cdf(dist, t); |
| |
| static const char* function = "madlib::modules::prob::cdf(" |
| "const students_t_distribution<%1%>&, %1%)"; |
| |
| RealType result; |
| if (!boost::math::detail::check_df(function, df, &result, Policy())) |
| return result; |
| |
| if (df >= 200) |
| return |
| boost::math::cdf( |
| boost::math::normal_distribution<RealType, Policy>(), |
| df >= 1000000 |
| ? t |
| : GleasonsNormalApproxForStudentsT(t, df) |
| ); |
| |
| // We first compute A = Pr[|T| < t] |
| RealType A = oneSidedStudentsT_CDF(t, static_cast<uint64_t>(df)); |
| |
| /* The Student-T distribution is obviously symmetric around t=0... */ |
| if (t < 0) |
| /* FIXME: If A is approximately 1, we will face a loss of significance. |
| * */ |
| return .5 * (1. - A); |
| else |
| /* While we only know A in [0,1] here, the end result will be in |
| * [0.5, 1]. Hence, there is no problem with adding 1 and A, even if |
| * A << 1. */ |
| return .5 * (1. + A); |
| } |
| |
| /** |
| * @brief Compute the complement of Student's cumulative distribution function |
| */ |
| template <class RealType, class Policy> |
| inline |
| RealType |
| cdf( |
| const boost::math::complemented2_type< |
| boost::math::students_t_distribution<RealType, Policy>, |
| RealType |
| >& c |
| ) { |
| RealType df = c.dist.degrees_of_freedom(); |
| if (df >= 200) { |
| static const char* function = "madlib::modules::prob::cdf(" |
| "const complement(students_t_distribution<%1%>&), %1%)"; |
| |
| RealType result; |
| if (!boost::math::detail::check_df(function, df, &result, Policy())) |
| return result; |
| |
| return |
| boost::math::cdf(complement( |
| boost::math::normal_distribution<RealType, Policy>(), |
| df >= 1000000 |
| ? c.param |
| : GleasonsNormalApproxForStudentsT(c.param, df) |
| )); |
| } |
| |
| return prob::cdf(c.dist, -c.param); |
| } |
| |
| template <class RealType, class Policy> |
| inline |
| RealType |
| pdf(const boost::math::students_t_distribution<RealType, Policy>& dist, |
| const RealType& t) { |
| return boost::math::pdf(dist, t); |
| } |
| |
| template <class RealType, class Policy> |
| inline |
| RealType |
| pdf( |
| const boost::math::complemented2_type< |
| boost::math::students_t_distribution<RealType, Policy>, |
| RealType |
| >& c |
| ) { |
| return boost::math::pdf(c); |
| } |
| |
| template <class RealType, class Policy> |
| inline |
| RealType |
| quantile(const boost::math::students_t_distribution<RealType, Policy>& dist, |
| const RealType& p) { |
| |
| using namespace boost::math; |
| |
| static const char* function = "madlib::modules::prob::quantile(" |
| "const students_t_distribution<%1%>&, %1%)"; |
| |
| // FIXME: Boost bug 6937 prevent proper argument validation. |
| // https://svn.boost.org/trac/boost/ticket/6937 |
| // Until this is fixed upstream, we do the following checks here. |
| RealType df = dist.degrees_of_freedom(); |
| RealType result; |
| if (!detail::check_df(function, df, &result, Policy()) |
| || !detail::check_probability(function, p, &result, Policy())) |
| return result; |
| |
| return boost::math::quantile(dist, p); |
| } |
| |
| template <class RealType, class Policy> |
| inline |
| RealType |
| quantile( |
| const boost::math::complemented2_type< |
| boost::math::students_t_distribution<RealType, Policy>, |
| RealType |
| >& c |
| ) { |
| return boost::math::quantile(c); |
| } |
| |
| } // namespace prob |
| |
| } // namespace modules |
| |
| } // namespace madlib |
| |
| #endif // defined(MADLIB_MODULES_PROB_STUDENT_T_HPP) |