/src/r-source/src/nmath/dt.c
Line | Count | Source |
1 | | /* |
2 | | * AUTHOR |
3 | | * Catherine Loader, catherine@research.bell-labs.com. |
4 | | * October 23, 2000. |
5 | | * |
6 | | * Merge in to R: |
7 | | * Copyright (C) 2000-2015 The R Core Team |
8 | | * |
9 | | * This program is free software; you can redistribute it and/or modify |
10 | | * it under the terms of the GNU General Public License as published by |
11 | | * the Free Software Foundation; either version 2 of the License, or |
12 | | * (at your option) any later version. |
13 | | * |
14 | | * This program is distributed in the hope that it will be useful, |
15 | | * but WITHOUT ANY WARRANTY; without even the implied warranty of |
16 | | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
17 | | * GNU General Public License for more details. |
18 | | * |
19 | | * You should have received a copy of the GNU General Public License |
20 | | * along with this program; if not, a copy is available at |
21 | | * https://www.R-project.org/Licenses/ |
22 | | * |
23 | | * |
24 | | * DESCRIPTION |
25 | | * |
26 | | * The t density is evaluated as |
27 | | * sqrt(n/2) / ((n+1)/2) * Gamma((n+3)/2) / Gamma((n+2)/2). |
28 | | * * (1+x^2/n)^(-n/2) |
29 | | * / sqrt( 2 pi (1+x^2/n) ) |
30 | | * |
31 | | * This form leads to a stable computation for all |
32 | | * values of n, including n -> 0 and n -> infinity. |
33 | | */ |
34 | | |
35 | | #include "nmath.h" |
36 | | #include "dpq.h" |
37 | | |
38 | | double dt(double x, double n, int give_log) |
39 | 0 | { |
40 | 0 | #ifdef IEEE_754 |
41 | 0 | if (ISNAN(x) || ISNAN(n)) |
42 | 0 | return x + n; |
43 | 0 | #endif |
44 | 0 | if (n <= 0) ML_WARN_return_NAN; |
45 | 0 | if(!R_FINITE(x)) |
46 | 0 | return R_D__0; |
47 | 0 | if(!R_FINITE(n)) |
48 | 0 | return dnorm(x, 0., 1., give_log); |
49 | | |
50 | 0 | double u, t = -bd0(n/2.,(n+1)/2.) + stirlerr((n+1)/2.) - stirlerr(n/2.), |
51 | 0 | x2n = x*x/n, // in [0, Inf] |
52 | 0 | ax = 0., // <- -Wpedantic |
53 | 0 | l_x2n; // := log(sqrt(1 + x2n)) = log(1 + x2n)/2 |
54 | 0 | bool lrg_x2n = (x2n > 1./DBL_EPSILON); |
55 | 0 | if (lrg_x2n) { // large x^2/n : |
56 | 0 | ax = fabs(x); |
57 | 0 | l_x2n = log(ax) - log(n)/2.; // = log(x2n)/2 = 1/2 * log(x^2 / n) |
58 | 0 | u = // log(1 + x2n) * n/2 = n * log(1 + x2n)/2 = |
59 | 0 | n * l_x2n; |
60 | 0 | } |
61 | 0 | else if (x2n > 0.2) { |
62 | 0 | l_x2n = log(1 + x2n)/2.; |
63 | 0 | u = n * l_x2n; |
64 | 0 | } else { |
65 | 0 | l_x2n = log1p(x2n)/2.; |
66 | 0 | u = -bd0(n/2.,(n+x*x)/2.) + x*x/2.; |
67 | 0 | } |
68 | | |
69 | | //old: return R_D_fexp(M_2PI*(1+x2n), t-u); |
70 | | |
71 | | // R_D_fexp(f,x) := (give_log ? -0.5*log(f)+(x) : exp(x)/sqrt(f)) |
72 | | // f = 2pi*(1+x2n) |
73 | | // ==> 0.5*log(f) = log(2pi)/2 + log(1+x2n)/2 = log(2pi)/2 + l_x2n |
74 | | // 1/sqrt(f) = 1/sqrt(2pi * (1+ x^2 / n)) |
75 | | // = 1/sqrt(2pi)/(|x|/sqrt(n)*sqrt(1+1/x2n)) |
76 | | // = M_1_SQRT_2PI * sqrt(n)/ (|x|*sqrt(1+1/x2n)) |
77 | 0 | if(give_log) |
78 | 0 | return t-u - (M_LN_SQRT_2PI + l_x2n); |
79 | | |
80 | | // else : if(lrg_x2n) : sqrt(1 + 1/x2n) ='= sqrt(1) = 1 |
81 | 0 | double I_sqrt_ = (lrg_x2n ? sqrt(n)/ax : exp(-l_x2n)); |
82 | 0 | return exp(t-u) * M_1_SQRT_2PI * I_sqrt_; |
83 | 0 | } |