stdtr.c
Go to the documentation of this file.
1 /* stdtr.c
2  *
3  * Student's t distribution
4  *
5  *
6  *
7  * SYNOPSIS:
8  *
9  * double t, stdtr();
10  * short k;
11  *
12  * y = stdtr( k, t );
13  *
14  *
15  * DESCRIPTION:
16  *
17  * Computes the integral from minus infinity to t of the Student
18  * t distribution with integer k > 0 degrees of freedom:
19  *
20  * t
21  * -
22  * | |
23  * - | 2 -(k+1)/2
24  * | ( (k+1)/2 ) | ( x )
25  * ---------------------- | ( 1 + --- ) dx
26  * - | ( k )
27  * sqrt( k pi ) | ( k/2 ) |
28  * | |
29  * -
30  * -inf.
31  *
32  * Relation to incomplete beta integral:
33  *
34  * 1 - stdtr(k,t) = 0.5 * incbet( k/2, 1/2, z )
35  * where
36  * z = k/(k + t**2).
37  *
38  * For t < -2, this is the method of computation. For higher t,
39  * a direct method is derived from integration by parts.
40  * Since the function is symmetric about t=0, the area under the
41  * right tail of the density is found by calling the function
42  * with -t instead of t.
43  *
44  * ACCURACY:
45  *
46  * Tested at random 1 <= k <= 25. The "domain" refers to t.
47  * Relative error:
48  * arithmetic domain # trials peak rms
49  * IEEE -100,-2 50000 5.9e-15 1.4e-15
50  * IEEE -2,100 500000 2.7e-15 4.9e-17
51  */
52 
53 /* stdtri.c
54  *
55  * Functional inverse of Student's t distribution
56  *
57  *
58  *
59  * SYNOPSIS:
60  *
61  * double p, t, stdtri();
62  * int k;
63  *
64  * t = stdtri( k, p );
65  *
66  *
67  * DESCRIPTION:
68  *
69  * Given probability p, finds the argument t such that stdtr(k,t)
70  * is equal to p.
71  *
72  * ACCURACY:
73  *
74  * Tested at random 1 <= k <= 100. The "domain" refers to p:
75  * Relative error:
76  * arithmetic domain # trials peak rms
77  * IEEE .001,.999 25000 5.7e-15 8.0e-16
78  * IEEE 10^-6,.001 25000 2.0e-12 2.9e-14
79  */
80 
81 
82 /*
83  * Cephes Math Library Release 2.3: March, 1995
84  * Copyright 1984, 1987, 1995 by Stephen L. Moshier
85  */
86 
87 #include "mconf.h"
88 #include <float.h>
89 
90 extern double MACHEP;
91 
92 double stdtr(int k, double t)
93 {
94  double x, rk, z, f, tz, p, xsqk;
95  int j;
96 
97  if (k <= 0) {
98  sf_error("stdtr", SF_ERROR_DOMAIN, NULL);
99  return (NAN);
100  }
101 
102  if (t == 0)
103  return (0.5);
104 
105  if (t < -2.0) {
106  rk = k;
107  z = rk / (rk + t * t);
108  p = 0.5 * incbet(0.5 * rk, 0.5, z);
109  return (p);
110  }
111 
112  /* compute integral from -t to + t */
113 
114  if (t < 0)
115  x = -t;
116  else
117  x = t;
118 
119  rk = k; /* degrees of freedom */
120  z = 1.0 + (x * x) / rk;
121 
122  /* test if k is odd or even */
123  if ((k & 1) != 0) {
124 
125  /* computation for odd k */
126 
127  xsqk = x / sqrt(rk);
128  p = atan(xsqk);
129  if (k > 1) {
130  f = 1.0;
131  tz = 1.0;
132  j = 3;
133  while ((j <= (k - 2)) && ((tz / f) > MACHEP)) {
134  tz *= (j - 1) / (z * j);
135  f += tz;
136  j += 2;
137  }
138  p += f * xsqk / z;
139  }
140  p *= 2.0 / M_PI;
141  }
142 
143 
144  else {
145 
146  /* computation for even k */
147 
148  f = 1.0;
149  tz = 1.0;
150  j = 2;
151 
152  while ((j <= (k - 2)) && ((tz / f) > MACHEP)) {
153  tz *= (j - 1) / (z * j);
154  f += tz;
155  j += 2;
156  }
157  p = f * x / sqrt(z * rk);
158  }
159 
160  /* common exit */
161 
162 
163  if (t < 0)
164  p = -p; /* note destruction of relative accuracy */
165 
166  p = 0.5 + 0.5 * p;
167  return (p);
168 }
169 
170 double stdtri(int k, double p)
171 {
172  double t, rk, z;
173  int rflg;
174 
175  if (k <= 0 || p <= 0.0 || p >= 1.0) {
176  sf_error("stdtri", SF_ERROR_DOMAIN, NULL);
177  return (NAN);
178  }
179 
180  rk = k;
181 
182  if (p > 0.25 && p < 0.75) {
183  if (p == 0.5)
184  return (0.0);
185  z = 1.0 - 2.0 * p;
186  z = incbi(0.5, 0.5 * rk, fabs(z));
187  t = sqrt(rk * z / (1.0 - z));
188  if (p < 0.5)
189  t = -t;
190  return (t);
191  }
192  rflg = -1;
193  if (p >= 0.5) {
194  p = 1.0 - p;
195  rflg = 1;
196  }
197  z = incbi(0.5 * rk, 0.5, 2.0 * p);
198 
199  if (DBL_MAX * z < rk)
200  return (rflg * INFINITY);
201  t = sqrt(rk / z - rk);
202  return (rflg * t);
203 }
stdtr
double stdtr(int k, double t)
Definition: stdtr.c:92
incbi
double incbi(double aa, double bb, double yy0)
Definition: incbi.c:51
stdtri
double stdtri(int k, double p)
Definition: stdtr.c:170
atan
const EIGEN_DEVICE_FUNC AtanReturnType atan() const
Definition: ArrayCwiseUnaryOps.h:283
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:119
j
std::ptrdiff_t j
Definition: tut_arithmetic_redux_minmax.cpp:2
incbet
double incbet(double aa, double bb, double xx)
Definition: incbet.c:283
MACHEP
double MACHEP
Definition: const.c:54
pybind_wrapper_test_script.z
z
Definition: pybind_wrapper_test_script.py:61
tree::f
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
Definition: testExpression.cpp:218
mconf.h
p
float * p
Definition: Tutorial_Map_using.cpp:9
sf_error
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Definition: sf_error.c:41
NULL
#define NULL
Definition: ccolamd.c:609
M_PI
#define M_PI
Definition: mconf.h:117
align_3::t
Point2 t(10, 10)
SF_ERROR_DOMAIN
@ SF_ERROR_DOMAIN
Definition: sf_error.h:16
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418


gtsam
Author(s):
autogenerated on Sat Nov 16 2024 04:05:00