From 8241956ddbf86f017acb3b4cc81d5d9334994720 Mon Sep 17 00:00:00 2001 From: Sean Pinkney Date: Tue, 11 Jul 2023 08:47:51 -0400 Subject: [PATCH] Update student_t_qf.stanfunctions (#85) * Update student_t_qf.stanfunctions * fix path * Formatting --------- Co-authored-by: GitHub Actions --- .../copula/student_t_copula.stanfunctions | 36 ++- .../distribution/student_t.stanfunctions | 41 +-- .../student_t_qf.stanfunctions | 263 +++++++++--------- 3 files changed, 172 insertions(+), 168 deletions(-) diff --git a/functions/copula/student_t_copula.stanfunctions b/functions/copula/student_t_copula.stanfunctions index ec7c586..6f85adb 100644 --- a/functions/copula/student_t_copula.stanfunctions +++ b/functions/copula/student_t_copula.stanfunctions @@ -1,4 +1,4 @@ - /** @addtogroup student_t_copula Student T Copula Functions +/** @addtogroup student_t_copula Student T Copula Functions * * @include \copula\student_t_copula.stanfunctions * @@ -15,23 +15,21 @@ * @param rho Real number [-1, 1] * @param nu Real \f$(0, +\infty)\f$ */ + +real bivariate_t_copula_lpdf(vector u, vector v, real rho, real nu) { + int N = num_elements(u); + vector[N] t1 = student_t_qf(u, nu); + vector[N] t2 = student_t_qf(v, nu); + vector[N] t1_sq = square(t1); + vector[N] t2_sq = square(t2); - real bivariate_t_copula_lpdf(vector u, vector v, real rho, real nu) { - int N = num_elements(u); - vector[N] t1 = student_t_qf(u, nu); - vector[N] t2 = student_t_qf(v, nu); - vector[N] t1_sq = square(t1); - vector[N] t2_sq = square(t2); - - real lpdf = N - * (-0.5 * log1m(rho ^ 2) + lgamma(0.5 * (nu + 2)) - + lgamma(0.5 * nu) - 2 * lgamma(0.5 * (nu + 1))); - lpdf += 0.5 * (nu + 1) - * sum(log1p(t1_sq / nu) + log1p(t2_sq / nu)); - lpdf += -0.5 * (nu + 2) - * sum(log1p((t1_sq - 2 * t1 .* t2 * rho + t2_sq) - / (nu * (1 - rho ^ 2)))); - return lpdf; - } + real lpdf = N + * (-0.5 * log1m(rho ^ 2) + lgamma(0.5 * (nu + 2)) + lgamma(0.5 * nu) + - 2 * lgamma(0.5 * (nu + 1))); + lpdf += 0.5 * (nu + 1) * sum(log1p(t1_sq / nu) + log1p(t2_sq / nu)); + lpdf += -0.5 * (nu + 2) + * sum(log1p((t1_sq - 2 * t1 .* t2 * rho + t2_sq) / (nu * (1 - rho ^ 2)))); + return lpdf; +} - /** @} */ \ No newline at end of file +/** @} */ \ No newline at end of file diff --git a/functions/distribution/student_t.stanfunctions b/functions/distribution/student_t.stanfunctions index 172550e..80e9f23 100644 --- a/functions/distribution/student_t.stanfunctions +++ b/functions/distribution/student_t.stanfunctions @@ -17,8 +17,8 @@ real student_t_lcdf_stan(real x, real df) { int lower_tail = 1; real lval; - - if (df <= 0) { + + if (df <= 0) { reject("df must be > 0. Found df = ", df); } @@ -26,25 +26,25 @@ real student_t_lcdf_stan(real x, real df) { return not_a_number(); } if (is_inf(df)) - return normal_lcdf(x | 0, 1); - + return normal_lcdf(x | 0, 1); + real nx = 1 + (x / df) * x; if (nx > 1e100) lval = -0.5 * df * (2 * log(abs(x)) - log(df)) - lbeta(0.5 * df, 0.5) - log(0.5 * df); - else - lval = df > x * x ? beta_lccdf(x * x / (df + x * x) | 0.5, df / 2) : beta_lcdf(1 / nx | df / 2, 0.5); + else + lval = df > x * x ? beta_lccdf(x * x / (df + x * x) | 0.5, df / 2) + : beta_lcdf(1 / nx | df / 2, 0.5); - if (x <= 0) { lower_tail = 0; } - - if (lower_tail == 1) { - return log1m(0.5 * exp(lval)); - } else { - return lval - log2(); - } - + + if (lower_tail == 1) { + return log1m(0.5 * exp(lval)); + } else { + return lval - log2(); + } + return lval; } @@ -59,7 +59,7 @@ real student_t_lcdf_stan(real x, real df) { */ real student_t_lccdf_stan(real x, real df) { real lval; - + if (df <= 0) { reject("df must be > 0. Found df = ", df); } @@ -69,14 +69,15 @@ real student_t_lccdf_stan(real x, real df) { } if (is_inf(df)) - return normal_lccdf(x | 0, 1); - + return normal_lccdf(x | 0, 1); + real nx = 1 + (x / df) * x; if (nx > 1e100) lval = -0.5 * df * (2 * log(abs(x)) - log(df)) - lbeta(0.5 * df, 0.5) - log(0.5 * df); - else - lval = df > x * x ? beta_lccdf(x * x / (df + x * x) | 0.5, df / 2) : beta_lcdf(1 / nx | df / 2, 0.5); - + else + lval = df > x * x ? beta_lccdf(x * x / (df + x * x) | 0.5, df / 2) + : beta_lcdf(1 / nx | df / 2, 0.5); + return x < 0 ? log1m(0.5 * exp(lval)) : lval - log2(); } diff --git a/functions/quantile_function/student_t_qf.stanfunctions b/functions/quantile_function/student_t_qf.stanfunctions index cda569f..0ef179c 100644 --- a/functions/quantile_function/student_t_qf.stanfunctions +++ b/functions/quantile_function/student_t_qf.stanfunctions @@ -1,4 +1,11 @@ -/** @ingroup qf +/** @addtogroup student_t_qf Student T Quantile functions + * + * @include \quantile_function\student_t_qf.stanfunctions + * + * \ingroup qf + * @{ */ + +/** * Student T Quantile Function * * @include \quantile_function\student_t_qf.stanfunctions @@ -9,7 +16,7 @@ * @return inverse CDF value * @throws reject if \f$ p \notin [0, 1] \f$ */ -#include univariate/student_t.stanfunctions +#include distribution/student_t.stanfunctions real student_t_qf(real p, real ndf) { if (is_nan(p) || p < 0 || p > 1) @@ -179,147 +186,144 @@ vector student_t_qf(vector u, real nu) { return x; } - /* Student T Log Quantile Function - * - * @include \quantile_function\student_t_qf.stanfunctions - * - * @author Sean Pinkney - * @param logp Real - * @param df Real \f$(0, +\infty)\f$ - * @return inverse CDF value - * @throws reject if \f$ p \notin [0, 1] \f$ - */ - real student_t_log_qf(real log_p, real ndf) { - real eps = 1e-12; - real d_epsilon = 2.220446e-16; - real d_max = 1.797693e+308; - real d_min = 2.225074e-308; - int d_mant_dig = 53; - real q; - real p = exp(log_p); - - if (is_nan(p) || is_nan(ndf)) { - return p + ndf; +/* Student T Log Quantile Function +* +* @include \quantile_function\student_t_qf.stanfunctions +* +* @author Sean Pinkney +* @param logp Real +* @param df Real \f$(0, +\infty)\f$ +* @return inverse CDF value +* @throws reject if \f$ p \notin [0, 1] \f$ +*/ +real student_t_log_qf(real log_p, real ndf) { + real eps = 1e-12; + real d_epsilon = 2.220446e-16; + real d_max = 1.797693e+308; + real d_min = 2.225074e-308; + int d_mant_dig = 53; + real q; + real p = exp(log_p); + + if (is_nan(p) || is_nan(ndf)) { + return p + ndf; + } + + if (ndf <= 0) { + reject("Invalid value for ndf"); + } + + if (ndf < 1) { + // Find the upper and lower bounds + real accu = 1e-13; + real Eps = 1e-11; + if (p > 1 - d_epsilon) { + return positive_infinity(); } - - if (ndf <= 0) { - reject("Invalid value for ndf"); + // real pp = min({1 - d_epsilon, p * (1 + Eps)}); + real log_pp = min({log1m(d_epsilon), log_p + log1p(Eps)}); + real ux = 1; + while (student_t_lcdf_stan(ux, ndf) < log_pp) { + ux *= 2; + } + log_pp = log_p - Eps; + real lx = -1; + while (student_t_lcdf_stan(lx, ndf) > log_pp) { + lx *= 2; } - if (ndf < 1) { - // Find the upper and lower bounds - real accu = 1e-13; - real Eps = 1e-11; - if (p > 1 - d_epsilon) { - return positive_infinity(); - } - // real pp = min({1 - d_epsilon, p * (1 + Eps)}); - real log_pp = min({log1m(d_epsilon), log_p + log1p(Eps)}); - real ux = 1; - while (student_t_lcdf_stan(ux , ndf) < log_pp) { - ux *= 2; - } - log_pp = log_p - Eps; - real lx = -1; - while (student_t_lcdf_stan(lx , ndf) > log_pp) { - lx *= 2; + // Find the quantile using interval halving + real nx = 0.5 * (lx + ux); + int iter = 0; + while ((ux - lx) / abs(nx) > accu && iter < 1000) { + iter = iter + 1; + if (student_t_lcdf_stan(nx, ndf) > log_p) { + ux = nx; + } else { + lx = nx; } - - // Find the quantile using interval halving - real nx = 0.5 * (lx + ux); - int iter = 0; - while ((ux - lx) / abs(nx) > accu && iter < 1000) { - iter = iter + 1; - if (student_t_lcdf_stan(nx , ndf) > log_p) { - ux = nx; - } else { - lx = nx; - } - nx = 0.5 * (lx + ux); + nx = 0.5 * (lx + ux); + } + return 0.5 * (lx + ux); + } + + if (ndf > 1e20) { + return std_normal_log_qf(log_p); + } + + int neg = log_p < -0.6931472 ? 1 : 0; + int is_neg_lower = neg; + real P = neg == 1 ? 2 * p : 2 * (0.5 - p + 0.5); + + P = min({max({P, 0}), 1}); + + if (abs(ndf - 2) < eps) { + if (P > d_min) { + if (3 * P < d_epsilon) { + q = 1 / sqrt(P); + } else if (P > 0.9) { + q = (1 - P) * sqrt(2 / (P * (2 - P))); + } else { + q = sqrt(2 / (P * (2 - P)) - 2); } - return 0.5 * (lx + ux); + } else { + q = positive_infinity(); } - - if (ndf > 1e20) { - return std_normal_log_qf(log_p); + } else if (ndf < 1. + eps) { + if (P == 1.) { + q = 0; + } else if (P > 0) { + q = 1 / tan(pi() * p / 2); + } else { + q = negative_infinity(); } + } else { + real x = 0; + real y = 0; + real log_P2 = 0; + real a = 1 / (ndf - 0.5); + real b = 48 / (a * a); + real c = ((20700 * a / b - 98) * a - 16) * a + 96.36; + real d = ((94.5 / (b + c) - 3) / b + 1) * sqrt(a * pi() / 2) * ndf; - int neg = log_p < -0.6931472 ? 1 : 0; - int is_neg_lower = neg; - real P = neg == 1 ? 2 * p : 2 * (0.5 - p + 0.5); - - P = min({max({P, 0}), 1}); - + y = pow((d * P), (2.0 / ndf)); + int P_ok = 0; + log_P2 = is_neg_lower == 1 ? log_p : log1m_exp(p); + x = (log(d) + log2() + log_P2) / ndf; + y = exp(2 * x); - if (abs(ndf - 2) < eps) { - if (P > d_min) { - if (3 * P < d_epsilon) { - q = 1 / sqrt(P); - } else if (P > 0.9) { - q = (1 - P) * sqrt(2 / (P * (2 - P))); - } else { - q = sqrt(2 / (P * (2 - P)) - 2); - } - } else { - q = positive_infinity(); - } - } else if (ndf < 1. + eps) { - if (P == 1.) { - q = 0; - } else if (P > 0) { - q = 1 / tan(pi() * p / 2); - } else { - q = negative_infinity(); + if ((ndf < 2.1 && P > 0.5) || y > 0.05 + a) { + x = inv_Phi(log_P2); + y = square(x); + + if (ndf < 5) { + c += 0.3 * (ndf - 4.5) * (x + 0.6); } + + c = (((0.05 * d * x - 5) * x - 7) * x - 2) * x + b + c; + y = (((((0.4 * y + 6.3) * y + 36) * y + 94.5) / c - y - 3) / b + 1) * x; + y = expm1(a * square(y)); + q = sqrt(ndf * y); + } else if (x < -log2() * d_mant_dig) { + q = sqrt(ndf) * exp(-x); } else { - real x = 0; - real y = 0; - real log_P2 = 0; - real a = 1 / (ndf - 0.5); - real b = 48 / (a * a); - real c = ((20700 * a / b - 98) * a - 16) * a + 96.36; - real d = ((94.5 / (b + c) - 3) / b + 1) * sqrt(a * pi() / 2) * ndf; + y = ((1 / (((ndf + 6) / (ndf * y) - 0.089 * d - 0.822) * (ndf + 2) * 3) + + 0.5 / (ndf + 4)) + * y - 1) + * (ndf + 1) / (ndf + 2) + 1 / y; - y = pow((d * P), (2.0 / ndf)); - int P_ok = 0; - log_P2 = is_neg_lower == 1 ? log_p : log1m_exp(p); - x = (log(d) + log2() + log_P2) / ndf; - y = exp(2 * x); - - if ((ndf < 2.1 && P > 0.5) || y > 0.05 + a) { - x = inv_Phi(log_P2); - y = square(x); - - if (ndf < 5) { - c += 0.3 * (ndf - 4.5) * (x + 0.6); - } - - c = (((0.05 * d * x - 5) * x - 7) * x - 2) * x + b + c; - y = (((((0.4 * y + 6.3) * y + 36) * y + 94.5) / c - y - 3) / b + 1) - * x; - y = expm1(a * square(y)); - q = sqrt(ndf * y); - } else if (x < -log2() * d_mant_dig) { - q = sqrt(ndf) * exp(-x); - } else { - y = ((1 - / (((ndf + 6) / (ndf * y) - 0.089 * d - 0.822) * (ndf + 2) * 3) - + 0.5 / (ndf + 4)) - * y - 1) - * (ndf + 1) / (ndf + 2) + 1 / y; - - q = sqrt(ndf * y); - } - } - - if (neg) { - q = -q; + q = sqrt(ndf * y); } - - return q; } + + if (neg) { + q = -q; + } + + return q; +} - vector student_t_log_qf(vector u, real nu) { +vector student_t_log_qf(vector u, real nu) { int N = num_elements(u); vector[N] x; @@ -328,4 +332,5 @@ vector student_t_qf(vector u, real nu) { } return x; -} \ No newline at end of file +} +/** @} */ \ No newline at end of file