SpecialFuncs.cpp
Go to the documentation of this file.
1 //==============================================================================
2 //
3 // This file is part of GNSSTk, the ARL:UT GNSS Toolkit.
4 //
5 // The GNSSTk is free software; you can redistribute it and/or modify
6 // it under the terms of the GNU Lesser General Public License as published
7 // by the Free Software Foundation; either version 3.0 of the License, or
8 // any later version.
9 //
10 // The GNSSTk is distributed in the hope that it will be useful,
11 // but WITHOUT ANY WARRANTY; without even the implied warranty of
12 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 // GNU Lesser General Public License for more details.
14 //
15 // You should have received a copy of the GNU Lesser General Public
16 // License along with GNSSTk; if not, write to the Free Software Foundation,
17 // Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110, USA
18 //
19 // This software was developed by Applied Research Laboratories at the
20 // University of Texas at Austin.
21 // Copyright 2004-2022, The Board of Regents of The University of Texas System
22 //
23 //==============================================================================
24 
25 //==============================================================================
26 //
27 // This software was developed by Applied Research Laboratories at the
28 // University of Texas at Austin, under contract to an agency or agencies
29 // within the U.S. Department of Defense. The U.S. Government retains all
30 // rights to use, duplicate, distribute, disclose, or release this software.
31 //
32 // Pursuant to DoD Directive 523024
33 //
34 // DISTRIBUTION STATEMENT A: This software has been approved for public
35 // release, distribution is unlimited.
36 //
37 //==============================================================================
38 
51 #include "SpecialFuncs.hpp"
52 #include <cmath>
53 #include <limits>
54 //#include "logstream.hpp" // TEMP
55 
56 namespace gnsstk
57 {
58 
59  /* Natural log of the gamma function for positive argument.
60  Gamma(x) = integral(0 to inf) { t^(x-1) exp(-t) dt }
61  param x argument, x must be > 0
62  return double ln(gamma(x)), the natural log of the gamma function of x.
63  throw if the input argument is <= 0 */
64  double lnGamma(double x)
65  {
66  try
67  {
68  static const double con[8] = {
69  76.18009172947146, -86.50532032941677, 24.01409824083091,
70  -1.231739572450155, 1.208650973866179e-3, -5.395239384953e-6,
71  1.000000000190015, 2.5066282746310005};
72 
73  if (x <= 0)
74  {
75  GNSSTK_THROW(Exception("Non-positive argument"));
76  }
77 
78  double y(x);
79  double t(x + 5.5);
80  t -= (x + 0.5) * ::log(t);
81  double s(con[6]);
82  for (int i = 0; i <= 5; i++)
83  s += con[i] / ++y;
84 
85  return (-t + ::log(con[7] * s / x));
86  }
87  catch (Exception& e)
88  {
89  GNSSTK_RETHROW(e);
90  }
91  }
92 
93  /* Gamma(x) the gamma function for positive argument.
94  Gamma(x) = integral(0 to inf) { t^(x-1) exp(-t) dt }
95  param x argument, x must be > 0
96  return double Gamma(x), the gamma function of x.
97  throw if the input argument is <= 0 */
98  double Gamma(double x)
99  {
100  try
101  {
102  return ::exp(lnGamma(double(x)));
103  }
104  catch (Exception& e)
105  {
106  GNSSTK_RETHROW(e);
107  }
108  }
109 
110  /* Factorial of an integer, returned as a double.
111  param n argument, n must be >= 0
112  return n! or factorial(n), as a double
113  throw if the input argument is < 0 */
114  double factorial(int n)
115  {
116  try
117  {
118  if (n < 0)
119  {
120  GNSSTK_THROW(Exception("Negative argument"));
121  }
122 
123  if (n > 32)
124  {
125  return ::exp(lnGamma(double(n + 1)));
126  }
127 
128  static double store[33] = {1.0, 1.0, 2.0, 6.0, 24.0, 120.0};
129  static int nstore = 5;
130 
131  while (nstore < n)
132  {
133  int i = nstore++;
134  store[nstore] = store[i] * nstore;
135  }
136 
137  return store[n];
138  }
139  catch (Exception& e)
140  {
141  GNSSTK_RETHROW(e);
142  }
143  }
144 
145  /* ln of Factorial of an integer, returned as a double.
146  param n argument, n must be >= 0
147  return ln(n!) or natural log of factorial(n), as a double
148  throw if the input argument is < 0 */
149  double lnFactorial(int n)
150  {
151  try
152  {
153  if (n < 0)
154  {
155  GNSSTK_THROW(Exception("Negative argument"));
156  }
157  if (n <= 1)
158  {
159  return 0.0;
160  }
161  return lnGamma(double(n + 1));
162  }
163  catch (Exception& e)
164  {
165  GNSSTK_RETHROW(e);
166  }
167  }
168 
169  /* Binomial coefficient (n k) = n!/[k!(n-k)!], 0 <= k <= n.
170  (n k) is the number of combinations of n things taken k at a time.
171  NB. (n+1 k) = [ (n+1)/(n-k+1) ] (n k) = (n k) + (n k-1)
172  NB. (n k+1) = [ (n-k)/(k+1) ] (n k)
173  param n int n must be >= 0
174  param k int k must be >= 0 and <= n
175  return (n k), the binomial coefficient
176  throw if the input argument do not satisfy 0 <= k <= n */
177  double binomialCoeff(int n, int k)
178  {
179  try
180  {
181  if (n < 0 || k > n)
182  {
183  GNSSTK_THROW(Exception("Invalid arguments"));
184  }
185 
186  if (n <= 32)
187  {
188  return (factorial(n) / (factorial(k) * factorial(n - k)));
189  }
190 
191  return floor(
192  0.5 + ::exp(lnFactorial(n) - lnFactorial(k) - lnFactorial(n - k)));
193  }
194  catch (Exception& e)
195  {
196  GNSSTK_RETHROW(e);
197  }
198  }
199 
200  /* Beta function. Beta(x,y)=Beta(y,x)=integral(0 to 1) {t^(x-1)*(1-t)^(y-1)
201  dt}. Also, Beta(x,y) = gamma(x)*gamma(y)/gamma(x+y). param x first
202  argument param y second argument return beta(x,y) throw if
203  either input argument is <= 0 */
204  double beta(double x, double y)
205  {
206  try
207  {
208  return ::exp(lnGamma(x) + lnGamma(y) - lnGamma(x + y));
209  }
210  catch (Exception& e)
211  {
212  GNSSTK_RETHROW(e);
213  }
214  }
215 
216  /* Incomplete gamma function P(a,x), evaluated using series representation.
217  P(a,x) = (1/gamma(a)) integral (0 to x) { exp(-t) t^(a-1) dt }
218  param a first argument, a > 0
219  param x second argument, x >= 0
220  return P(a,x)
221  throw if input arguments have a <= 0 or x < 0 */
222  double seriesIncompGamma(double a, double x)
223  {
224  try
225  {
226  if (x < 0)
227  {
228  GNSSTK_THROW(Exception("Negative first argument"));
229  }
230  if (a <= 0)
231  {
232  GNSSTK_THROW(Exception("Non-positive second argument"));
233  }
234 
235  static const int imax(1000);
236  static const double eps(10 * std::numeric_limits<double>().epsilon());
237 
238  double lngamma(lnGamma(a));
239 
240  double atmp(a), sum(1.0 / a);
241  double del(sum);
242  for (int i = 1; i <= imax; i++)
243  {
244  ++atmp;
245  del *= x / atmp;
246  sum += del;
247  if (::fabs(del) < ::fabs(sum) * eps)
248  {
249  return (sum * ::exp(-x + a * ::log(x) - lngamma));
250  }
251  }
252  GNSSTK_THROW(Exception("Overflow; first arg too big"));
253  }
254  catch (Exception& e)
255  {
256  GNSSTK_RETHROW(e);
257  }
258 
259  return 0.0;
260  }
261 
262  /* Incomplete gamma function Q(a,x), evaluated using continued fractions.
263  Q(a,x) = (1/gamma(a)) integral (x to inf) { exp(-t) t^(a-1) dt }
264  param a first argument, a > 0
265  param x second argument, x >= 0
266  return Q(a,x)
267  throw if input arguments have a <= 0 or x < 0 */
268  double contfracIncompGamma(double a, double x)
269  {
270  try
271  {
272  if (x < 0)
273  {
274  GNSSTK_THROW(Exception("Negative first argument"));
275  }
276  if (a <= 0)
277  {
278  GNSSTK_THROW(Exception("Non-positive second argument"));
279  }
280 
281  static const int imax(1000);
282  static const double eps(10 * std::numeric_limits<double>().epsilon());
283  static const double fpmin(10 * std::numeric_limits<double>().min());
284 
285  double lngamma(lnGamma(a));
286 
287  double b(x + 1.0 - a), c(1.0 / fpmin);
288  double d(1.0 / b);
289  double h(d);
290  int i;
291  for (i = 1; i <= imax; i++)
292  {
293  double an(-i * (i - a));
294  b += 2.0;
295  d = an * d + b;
296  if (::fabs(d) < fpmin)
297  {
298  d = fpmin;
299  }
300  c = b + an / c;
301  if (::fabs(c) < fpmin)
302  {
303  c = fpmin;
304  }
305  d = 1.0 / d;
306  double del(d * c);
307  h *= del;
308  if (::fabs(del - 1.0) < eps)
309  {
310  break;
311  }
312  }
313 
314  if (i > imax)
315  {
316  GNSSTK_THROW(Exception("Overflow; first arg too big"));
317  }
318 
319  return (::exp(-x + a * ::log(x) - lngamma) * h);
320  }
321  catch (Exception& e)
322  {
323  GNSSTK_RETHROW(e);
324  }
325  }
326 
327  /* Incomplete gamma function P(a,x), a,x > 0.
328  P(a,x) = (1/gamma(a)) integral (0 to x) { exp(-t) t^(a-1) dt }; a > 0, x
329  >= 0 param a first argument, a > 0 param x second argument, x >= 0
330  return P(a,x)
331  throw if input arguments have a <= 0 or x < 0 */
332  double incompGamma(double a, double x)
333  {
334  try
335  {
336  if (x < 0)
337  {
338  GNSSTK_THROW(Exception("Negative first argument"));
339  }
340  if (a <= 0)
341  {
342  GNSSTK_THROW(Exception("Non-positive second argument"));
343  }
344 
345  if (x < a + 1.0)
346  {
347  return seriesIncompGamma(a, x);
348  }
349  else
350  {
351  return (1.0 - contfracIncompGamma(a, x));
352  }
353  }
354  catch (Exception& e)
355  {
356  GNSSTK_RETHROW(e);
357  }
358  }
359 
360  /* Complement of incomplete gamma function Q(a,x), a > 0, x >= 0.
361  Q(a,x) = (1/gamma(a)) integral (x to inf) { exp(-t) t^(a-1) dt }
362  param a first argument, a > 0
363  param x second argument, x >= 0
364  return Q(a,x)
365  throw if input arguments have a <= 0 or x < 0 */
366  double compIncompGamma(double a, double x)
367  {
368  try
369  {
370  if (x < 0)
371  {
372  GNSSTK_THROW(Exception("Negative first argument"));
373  }
374  if (a <= 0)
375  {
376  GNSSTK_THROW(Exception("Non-positive second argument"));
377  }
378 
379  if (x < a + 1.0)
380  {
381  return (1.0 - seriesIncompGamma(a, x));
382  }
383  else
384  {
385  return contfracIncompGamma(a, x);
386  }
387  }
388  catch (Exception& e)
389  {
390  GNSSTK_RETHROW(e);
391  }
392  }
393 
394  /* Error function erf(x). erf(x) = 2/sqrt(pi)*integral (0 to x) {exp(-t^2) dt}
395  param x input argument return erf(x) */
396  double errorFunc(double x)
397  {
398  if (x < 0)
399  {
400  GNSSTK_THROW(Exception("Negative first argument"));
401  }
402  try
403  {
404  return (x < 0.0 ? -incompGamma(0.5, x * x) : incompGamma(0.5, x * x));
405  }
406  catch (Exception& e)
407  {
408  GNSSTK_RETHROW(e);
409  }
410  }
411 
412  /* Complementary error function erfc(x). erfc(x) = 1-erf(x)
413  param x input argument
414  return erfc(x) */
415  double compErrorFunc(double x)
416  {
417  if (x < 0)
418  {
419  GNSSTK_THROW(Exception("Negative first argument"));
420  }
421  try
422  {
423  return (x < 0.0 ? 1.0 + incompGamma(0.5, x * x)
424  : compIncompGamma(0.5, x * x));
425  }
426  catch (Exception& e)
427  {
428  GNSSTK_RETHROW(e);
429  }
430  }
431 
432  /* Compute continued fractions portion of incomplete beta function I_x(a,b)
433  Routine used internally for Incomplete beta function I_x(a,b) */
434  double cfIBeta(double x, double a, double b)
435  {
436  static const int imax(1000);
437  static const double eps(10 * std::numeric_limits<double>().epsilon());
438  static const double fpmin(10 * std::numeric_limits<double>().min());
439  const double qab(a + b);
440  const double qap(a + 1.0);
441  const double qam(a - 1.0);
442  double c(1), d(1 - qab * x / qap), aa, del;
443  if (::fabs(d) < fpmin)
444  {
445  d = fpmin;
446  }
447  d = 1.0 / d;
448  double h(d);
449  int i, i2;
450  for (i = 1; i <= imax; i++)
451  {
452  i2 = 2 * i;
453  aa = i * (b - i) * x / ((qam + i2) * (a + i2));
454  d = 1.0 + aa * d;
455  if (::fabs(d) < fpmin)
456  {
457  d = fpmin;
458  }
459  c = 1.0 + aa / c;
460  if (::fabs(c) < fpmin)
461  {
462  c = fpmin;
463  }
464  d = 1.0 / d;
465  h *= d * c;
466  aa = -(a + i) * (qab + i) * x / ((a + i2) * (qap + i2));
467  d = 1.0 + aa * d;
468  if (::fabs(d) < fpmin)
469  {
470  d = fpmin;
471  }
472  c = 1.0 + aa / c;
473  if (::fabs(c) < fpmin)
474  {
475  c = fpmin;
476  }
477  d = 1.0 / d;
478  del = d * c;
479  h *= del;
480  if (::fabs(del - 1.0) < eps)
481  {
482  break;
483  }
484  }
485  if (i > imax)
486  {
487  GNSSTK_THROW(Exception("Overflow; a or b too big"));
488  }
489 
490  return h;
491  }
492 
493  /* Incomplete beta function I_x(a,b), 0<=x<=1, a,b>0
494  I sub x (a,b) = (1/beta(a,b)) integral (0 to x) { t^(a-1)*(1-t)^(b-1)dt }
495  param x input value, 0 <= x <= 1
496  param a input value, a > 0
497  param b input value, b > 0
498  return Incomplete beta function I_x(a,b) */
499  double incompleteBeta(double x, double a, double b)
500  {
501  if (x < 0 || x > 1)
502  {
503  GNSSTK_THROW(Exception("Invalid x argument"));
504  }
505  if (a <= 0 || b <= 0)
506  {
507  GNSSTK_THROW(Exception("Non-positive argument"));
508  }
509 
510  if (x == 0)
511  {
512  return 0.0;
513  }
514  if (x == 1)
515  {
516  return 1.0;
517  }
518 
519  try
520  {
521  double factor = ::exp(lnGamma(a + b) - lnGamma(a) - lnGamma(b) +
522  a * ::log(x) + b * ::log(1.0 - x));
523  if (x < (a + 1.0) / (a + b + 2.0))
524  {
525  return factor * cfIBeta(x, a, b) / a;
526  }
527  else
528  {
529  return 1.0 - factor * cfIBeta(1.0 - x, b, a) / b;
530  }
531  }
532  catch (Exception& e)
533  {
534  GNSSTK_RETHROW(e);
535  }
536  }
537 
538  // ----------------- probability distributions -----------------------
539 
540  /* Normal distribution of sample mean mu and sample std deviation sig
541  (location and scale parameters, resp.).
542  \code
543  NormalPDF(x,mu,sig) = exp(-(x-mu)*(x-mu)/(2*sig*sig));
544  NormalCDF(x,mu,sig) = 0.5*(1+erf((x-mu)/(::sqrt(2)*sig));
545  \endcode
546  For both theoretical and practical reasons, the normal distribution is
547  probably the most important distribution in statistics.
548  Many classical statistical tests are based on the assumption that the data
549  follow a normal distribution. (This assumption should be tested before
550  applying these tests.) In modeling applications, such as linear and
551  non-linear regression, the error term is often assumed to follow a normal
552  distribution with fixed location (mu) and scale (sig). The normal
553  distribution is widely used. Part of its appeal is that it is well behaved
554  and mathematically tractable. However, the central limit theorem provides
555  a theoretical basis for why it has wide applicability. The central limit
556  theorem states that as the sample size n becomes large, the following
557  occur:
558  The sampling distribution of the mean becomes approximately normal
559  regardless of the distribution of the original variable.
560  The sampling distribution of the mean is centered at the population
561  mean,
562  mu, of the original variable. In addition, the standard deviation of
563  the sampling distribution of the mean approaches sig/sqrt(n).
564  Probability density function (PDF) of the Normal distribution.
565  Ref http://www.itl.nist.gov/div898/handbook/ 1.3.6.6.1
566  param x input statistic
567  param mu mean of the sample (location parameter of the distribution)
568  param sig std dev of the sample (scale parameter of the distribution)
569  return Normal distribution probability density */
570  double NormalPDF(double x, double mu, double sig)
571  {
572  try
573  {
574  return (::exp(-(x - mu) * (x - mu) / (2.0 * sig * sig)));
575  }
576  catch (Exception& e)
577  {
578  GNSSTK_RETHROW(e);
579  }
580  }
581 
582  /* Cumulative distribution function (CDF) of the Normal-distribution.
583  Ref http://www.itl.nist.gov/div898/handbook/ 1.3.6.6.1
584  param x input statistic
585  param mu mean of the sample (location parameter of the distribution)
586  param sig std dev of the sample (scale parameter of the distribution)
587  return Normal distribution probability */
588  double NormalCDF(double x, double mu, double sig)
589  {
590  if (sig <= 0.0)
591  {
592  GNSSTK_THROW(Exception("Non-positive sigma"));
593  }
594 
595  try
596  {
597  static const double sqrt2(::sqrt(2.0));
598  double arg(x - mu);
599  double erf = errorFunc(::fabs(arg) / (sqrt2 * sig));
600  return (0.5 * (1.0 + (arg < 0.0 ? -erf : erf)));
601  }
602  catch (Exception& e)
603  {
604  GNSSTK_RETHROW(e);
605  }
606  }
607 
608  /* Normal-distribution percent point function, or inverse of the Normal CDF.
609  This function(prob,mu,sig) == X where prob = NormalCDF(X,mu,sig).
610  Ref http://www.itl.nist.gov/div898/handbook/ 1.3.6.6.1
611  param prob probability or significance level of the test, >=0 and < 1
612  param mu mean of the sample (location parameter of the distribution)
613  param sig std dev of the sample (scale parameter of the distribution)
614  return X the statistic at this probability */
615  double invNormalCDF(double prob, double mu, double sig)
616  {
617  try
618  {
619  if (prob < 0 || prob >= 1)
620  {
621  GNSSTK_THROW(Exception("Invalid probability argument"));
622  }
623  if (sig <= 0.0)
624  {
625  GNSSTK_THROW(Exception("Non-positive sigma"));
626  }
627 
628  static const double eps(1000000 *
629  std::numeric_limits<double>().epsilon());
630 
631  if (prob < eps)
632  {
633  return 0.0;
634  }
635  if (1.0 - prob < eps)
636  {
637  GNSSTK_THROW(Exception("Invalid probability -- too close to 1.0"));
638  }
639 
640  /* find X such that prob == NormalCDF(X,mu,sig); use bracket method
641  we know 0.5 = NormalCDF(muwhen prob = 0.5, X = mu
642  also invNormalCDF(1-prob,mu,sig) = 2*mu - invNormalCDF(prob,mu,sig)
643  so make alpha >= 0.5 */
644  bool swap(false);
645  double alpha(prob);
646  if (prob < 0.5)
647  {
648  swap = true;
649  alpha = 1.0 - prob;
650  }
651 
652  // we know a0 = NormalCDF(X0,mu,sig) where a0=0.5,X0=mu and alpha >= 0.5
653  double X, X0(mu), X1, a;
654  // first find X1 such that a1 = NormalCDF(X1,mu,sig) and a1 > alpha
655  X1 = 2.0;
656  while ((a = NormalCDF(X1, mu, sig)) <= alpha)
657  {
658  X1 *= 2.0;
659  }
660 
661  // bracket
662  int niter(0); // iteration count
663  while (1)
664  {
665  X = (X0 + X1) / 2.0;
666  a = NormalCDF(X, mu, sig);
667  // LOG(INFO) << "LOOP invNormalCDF X = " << niter << " " <<
668  // std::fixed
669  //<< std::setprecision(15) << X0 << " < " << X << " < " << X1
670  //<< " and a = " << a << " alpha = " << alpha << " a-alpha = "
671  //<< std::scientific << std::setprecision(2) << a-alpha << " eps "
672  //<< eps
673  //<< " fabs(a-alpha)-eps " << ::fabs(alpha-a)-eps;
674  if (::fabs(alpha - a) < eps)
675  {
676  break;
677  }
678  if (a > alpha)
679  {
680  X1 = X;
681  }
682  else
683  {
684  X0 = X;
685  }
686  if (++niter > 1000)
687  {
688  GNSSTK_THROW(Exception("Failed to converge"));
689  }
690  }
691 
692  return (swap ? 2.0 * mu - X : X);
693  }
694  catch (Exception& e)
695  {
696  GNSSTK_RETHROW(e);
697  }
698  }
699 
700  /* Probability density function (PDF) of the Chi-square distribution.
701  The chi-square distribution results when n independent variables with
702  standard normal distributions are squared and summed; x=RSS(variables).
703 
704  A chi-square test (Snedecor and Cochran, 1983) can be used to test if the
705  standard deviation of a population is equal to a specified value. This
706  test can be either a two-sided test or a one-sided test. The two-sided
707  version tests against the alternative that the true standard deviation is
708  either less than or greater than the specified value. The one-sided
709  version only tests in one direction. The chi-square hypothesis test is
710  defined as:
711  H0: sigma = sigma0
712  Ha: sigma < sigma0 for a lower one-tailed test
713  sigma > sigma0 for an upper one-tailed test
714  sigma <>sigma0 for a two-tailed test
715  Test Statistic: T = T = (N-1)*(s/sigma0)**2
716  where N is the sample size and s is the sample standard deviation.
717  The key element of this formula is the ratio s/sigma0 which compares the
718  ratio of the sample standard deviation to the target standard deviation.
719  As this ratio deviates from 1, the more likely is rejection of the null
720  hypothesis. Significance Level: alpha. Critical Region: Reject the null
721  hypothesis that the standard deviation is a specified value, sigma0, if
722  T > chisquare(alpha,N-1) for an upper one-tailed alternative
723  T < chisquare(1-alpha,N-1) for a lower one-tailed alternative
724  T < chisquare(1-alpha,N-1) for a two-tailed test or
725  T < chisquare(1-alpha,N-1)
726  where chi-square(p,N-1) is the critical value or inverseCDF of the
727  chi-square distribution with N-1 degrees of freedom.
728 
729  param x input statistic, equal to an RSS(); x >= 0
730  param n input value for number of degrees of freedom, n > 0
731  return probability Chi-square probability (xsq,n) */
732  double ChisqPDF(double x, int n)
733  {
734  if (x < 0)
735  {
736  GNSSTK_THROW(Exception("Negative statistic"));
737  }
738  if (n <= 0)
739  {
740  GNSSTK_THROW(Exception("Non-positive degrees of freedom"));
741  }
742 
743  try
744  {
745  double dn(double(n) / 2.0);
746  return (::exp(-x / 2.0) * ::pow(x, dn - 1.0) /
747  (::pow(2.0, dn) * Gamma(dn)));
748  }
749  catch (Exception& e)
750  {
751  GNSSTK_RETHROW(e);
752  }
753  }
754 
755  /* Cumulative distribution function (CDF) of the Chi-square-distribution.
756  Ref http://www.itl.nist.gov/div898/handbook/ 1.3.6.6.6
757  param x input statistic value, the RSS of variances, X >= 0
758  param n degrees of freedom of sample, n > 0
759  return probability that the sample variance is less than X. */
760  double ChisqCDF(double x, int n)
761  {
762  if (x < 0)
763  {
764  GNSSTK_THROW(Exception("Negative statistic"));
765  }
766  if (n <= 0)
767  {
768  GNSSTK_THROW(Exception("Non-positive degrees of freedom"));
769  }
770 
771  try
772  {
773  /* NB this incompGamma(n/2,x/2) == NIST's
774  incompGamma(n/2,x/2)/Gamma(n/2) */
775  return incompGamma(double(n) / 2.0, x / 2.0);
776  }
777  catch (Exception& e)
778  {
779  GNSSTK_RETHROW(e);
780  }
781  }
782 
783  /* Chi-square-distribution percent point function, or inverse of the Chisq
784  CDF. This function(alpha,N) == Y where alpha = ChisqCDF(Y,N). Ref
785  http://www.itl.nist.gov/div898/handbook/ 1.3.6.6.6 param alpha probability
786  or significance level of the test, >=0 and < 1 param n degrees of
787  freedom of sample, n > 0 return X the statistic (an RSS of
788  variances) at this probability */
789  double invChisqCDF(double alpha, int n)
790  {
791  try
792  {
793  if (alpha < 0 || alpha >= 1)
794  {
795  GNSSTK_THROW(Exception("Invalid probability argument"));
796  }
797  if (n <= 0)
798  {
799  GNSSTK_THROW(Exception("Non-positive degree of freedom"));
800  }
801 
802  static const double eps(1000000 *
803  std::numeric_limits<double>().epsilon());
804  if (alpha < eps)
805  {
806  return 0.0;
807  }
808  if (1.0 - alpha < eps)
809  {
810  GNSSTK_THROW(Exception("Invalid probability -- too close to 1.0"));
811  }
812 
813  /* find X such that alpha == ChisqCDF(X,n); use bracket method
814  we know a0 = ChisqCDF(F0,n) where a0=F0=0 */
815  double X, X0(0.0), X1, a;
816  // first find X1 such that a1 = ChisqCDF(X1,N) and a1 > alpha
817  X1 = 2.0;
818  while ((a = ChisqCDF(X1, n)) <= alpha)
819  {
820  X1 *= 2.0;
821  }
822 
823  // bracket
824  int niter(0); // iteration count
825  while (1)
826  {
827  X = (X0 + X1) / 2.0;
828  a = ChisqCDF(X, n);
829  // LOG(INFO) << "LOOP invChisqCDF X = " << niter << " " <<
830  // std::fixed
831  //<< std::setprecision(15) << X0 << " < " << X << " < " << X1
832  //<< " and a = " << a << " alpha = " << alpha << " a-alpha = "
833  //<< std::scientific << std::setprecision(2) << a-alpha << " eps "
834  //<< eps
835  //<< " fabs(a-alpha)-eps " << ::fabs(alpha-a)-eps;
836  if (::fabs(alpha - a) < eps)
837  {
838  break;
839  }
840  if (a > alpha)
841  {
842  X1 = X;
843  }
844  else
845  {
846  X0 = X;
847  }
848  if (++niter > 1000)
849  {
850  GNSSTK_THROW(Exception("Failed to converge"));
851  }
852  }
853 
854  return X;
855  }
856  catch (Exception& e)
857  {
858  GNSSTK_RETHROW(e);
859  }
860  }
861 
862  /* Probability density function (PDF) of the Student t distribution.
863  The null hypotheses that test the true mean, mu, against the standard or
864  assumed mean, mu0 are:
865  H0: mu = mu0
866  H0: mu <= mu0
867  H0: mu >= mu0
868  The basic statistics for the test are the sample mean and the standard
869  deviation. The form of the test statistic depends on whether the poulation
870  standard deviation, sigma, is known or is estimated from the data at hand.
871  The more typical case is where the standard deviation must be estimated
872  from the data, and the test statistic is
873  t = (Ybar - mu0/(s/SQRT(N))
874  where the sample mean is
875  Ybar = (1/N)*SUM[i=1 to N]Y(i)
876  and the sample standard deviation is
877  s = SQRT{(1/(N-1))*SUM[i=1 to N][Y(i) - Ybar)**2}
878  with N - 1 degrees of freedom.
879  For a test at significance level (probability) alpha, where alpha is
880  chosen to be small, typically .01, .05 or .10, the hypothesis associated
881  with each case enumerated above is rejected if:
882  |t| >= t(alpha/2,N-1)
883  t >= t(alpha,N-1)
884  t <= -t(alpha,N-1)
885  where t(alpha/2,N-1) is the upper alpha/2 critical value (inverse CDF)
886  from the t distribution with N-1 degrees of freedom.
887  param X input statistic
888  param n input value for number of degrees of freedom, n > 0
889  return probability density */
890  double StudentsPDF(double X, int n)
891  {
892  if (n <= 0)
893  {
894  GNSSTK_THROW(Exception("Non-positive degrees of freedom"));
895  }
896 
897  try
898  {
899  double dn(n);
900  return (::pow(1.0 + X * X / dn, -(dn + 1) / 2.0) /
901  (::sqrt(dn) * beta(0.5, 0.5 * dn)));
902  }
903  catch (Exception& e)
904  {
905  GNSSTK_RETHROW(e);
906  }
907  }
908 
909  /* Cumulative Distribution Function CDF() for Student-t-distribution CDF.
910  If X is a random variable following a normal distribution with mean zero
911  and variance unity, and chisq is a random variable following an
912  independent chi-square distribution with n degrees of freedom, then the
913  distribution of the ratio X/sqrt(chisq/n) is called Student's
914  t-distribution with n degrees of freedom. The probability that
915  |X/sqrt(chisq/n)| will be less than a fixed constant t is StudentCDF(t,n);
916  Ref http://www.itl.nist.gov/div898/handbook/ 1.3.6.6.4
917  Abramowitz and Stegun 26.7.1
918  param t input statistic value
919  param n degrees of freedom of first sample, n > 0
920  return probability that the sample is less than X. */
921  double StudentsCDF(double t, int n)
922  {
923  if (n <= 0)
924  {
925  GNSSTK_THROW(Exception("Non-positive degree of freedom"));
926  }
927 
928  try
929  {
930  // NB StudentsCDF(-t,n) = 1.0-StudentsCDF(t,n);
931  double x = 0.5 * incompleteBeta(double(n) / (t * t + double(n)),
932  double(n) / 2, 0.5);
933  if (t >= 0.0)
934  {
935  return (1.0 - x);
936  }
937  return (x);
938  }
939  catch (Exception& e)
940  {
941  GNSSTK_RETHROW(e);
942  }
943  }
944 
945  /* Students-t-distribution percent point function, or inverse of the Student
946  CDF. This function(prob,n) == Y where prob = StudentsCDF(Y,n). Ref
947  http://www.itl.nist.gov/div898/handbook/ 1.3.6.6.4 param prob probability
948  or significance level of the test, >=0 and < 1 param n degrees of
949  freedom of sample, n > 0 return t the statistic at this probability */
950  double invStudentsCDF(double prob, int n)
951  {
952  try
953  {
954  if (prob < 0 || prob >= 1)
955  {
956  GNSSTK_THROW(Exception("Invalid probability argument"));
957  }
958  if (n <= 0)
959  {
960  GNSSTK_THROW(Exception("Non-positive degree of freedom"));
961  }
962 
963  static const double eps(1000000 *
964  std::numeric_limits<double>().epsilon());
965  if (prob < eps)
966  {
967  return 0.0;
968  }
969  if (1.0 - prob < eps)
970  {
971  GNSSTK_THROW(Exception("Invalid probability -- too close to 1.0"));
972  }
973 
974  // find X such that prob == StudentsCDF(X,n); use bracket method
975 
976  // NB StudentsCDF(-t,n) = 1.0-StudentsCDF(t,n);
977  bool swap(false);
978  double alpha(prob);
979  if (prob < 0.5)
980  {
981  swap = true;
982  alpha = 1.0 - prob;
983  }
984 
985  // we know a0 = StudentsCDF(t0,n) where a0=0.5,t0=0 and alpha >= 0.5
986  double t, t0(0.0), t1, a;
987  // first find t1 such that a1 = StudentsCDF(t1,n) and a1 > alpha
988  t1 = 2.0;
989  while ((a = StudentsCDF(t1, n)) <= alpha)
990  {
991  t1 *= 2.0;
992  }
993 
994  // bracket
995  int niter(0); // iteration count
996  while (1)
997  {
998  t = (t0 + t1) / 2.0;
999  a = StudentsCDF(t, n);
1000  // LOG(INFO) << "LOOP invStudentsCDF t = " << niter << " " <<
1001  // std::fixed
1002  //<< std::setprecision(15) << t0 << " < " << t << " < " << t1
1003  //<< " and a = " << a << " alpha = " << alpha << " a-alpha = "
1004  //<< std::scientific << std::setprecision(2) << a-alpha << " eps "
1005  //<< eps
1006  //<< " fabs(a-alpha)-eps " << ::fabs(alpha-a)-eps;
1007  if (::fabs(alpha - a) < eps)
1008  {
1009  break;
1010  }
1011  if (a > alpha)
1012  {
1013  t1 = t;
1014  }
1015  else
1016  {
1017  t0 = t;
1018  }
1019  if (++niter > 1000)
1020  {
1021  GNSSTK_THROW(Exception("Failed to converge"));
1022  }
1023  }
1024 
1025  return (swap ? -t : t);
1026  }
1027  catch (Exception& e)
1028  {
1029  GNSSTK_RETHROW(e);
1030  }
1031  }
1032 
1033  /* F-distribution cumulative distribution function FDistCDF(F,n1,n2) F>=0
1034  n1,n2>0. This function occurs in the statistical test of whether two
1035  observed samples have the same variance. If F is the ratio of the observed
1036  dispersion (variance) of the first sample to that of the second, where the
1037  first sample has n1 degrees of freedom and the second has n2 degrees of
1038  freedom, then this function returns the probability that F would be as
1039  large as it is if the first sample's distribution has smaller variance
1040  than the second's. In other words, FDistCDF(f,n1,n2) is the significance
1041  level at which the hypothesis "sample 1 has smaller variance than sample
1042  2" can be rejected. A small numerical value implies a significant
1043  rejection, in turn implying high confidence in the hypothesis "sample 1
1044  has variance greater than or equal to that of sample 2". Ref
1045  http://www.itl.nist.gov/div898/handbook/ 1.3.6.6.5 param F input
1046  statistic value, the ratio variance1/variance2, F >= 0 param n1 degrees
1047  of freedom of first sample, n1 > 0 param n2 degrees of freedom of
1048  second sample, n2 > 0 return probability that the sample is less
1049  than F. */
1050  double FDistCDF(double F, int n1, int n2)
1051  {
1052  if (F < 0)
1053  {
1054  GNSSTK_THROW(Exception("Negative statistic"));
1055  }
1056  if (n1 <= 0 || n2 <= 0)
1057  {
1058  GNSSTK_THROW(Exception("Non-positive degree of freedom"));
1059  }
1060 
1061  try
1062  {
1063  return (1.0 -
1064  incompleteBeta(double(n2) / (double(n2) + double(n1) * F),
1065  double(n2) / 2.0, double(n1) / 2.0));
1066  }
1067  catch (Exception& e)
1068  {
1069  GNSSTK_RETHROW(e);
1070  }
1071  }
1072 
1073  /* Probability density function for F distribution
1074  The F distribution is the ratio of two chi-square distributions with
1075  degrees of freedom N1 and N2, respectively, where each chi-square has
1076  first been divided by its degrees of freedom. An F-test (Snedecor and
1077  Cochran, 1983) is used to test if the standard deviations of two
1078  populations are equal. This test can be a two-tailed test or a one-tailed
1079  test. The F hypothesis test is defined as:
1080  H0: s1 = s2 (sN is sigma or std deviation)
1081  Ha: s1 < s2 for a lower one tailed test
1082  s1 > s2 for an upper one tailed test
1083  s1 != s2 for a two tailed test
1084  Test Statistic: F = s1^2/s2^2 where s1^2 and s2^2 are the sample
1085  variances. The more this ratio deviates from 1, the stronger the evidence
1086  for unequal population variances. Significance Level is alpha, a
1087  probability (0<=alpha<=1). The hypothesis that the two standard deviations
1088  are equal is rejected if
1089  F > PP(alpha,N1-1,N2-1) for an upper one-tailed test
1090  F < PP(1-alpha,N1-1,N2-1) for a lower one-tailed test
1091  F < PP(1-alpha/2,N1-1,N2-1) for a two-tailed test
1092  F > PP(alpha/2,N1-1,N2-1)
1093  where PP(alpha,k-1,N-1) is the percent point function of the F
1094  distribution [PPfunc is inverse of the CDF : PP(alpha,N1,N2) == F where
1095  alpha=CDF(F,N1,N2)] with N1 and N2 degrees of freedom and a significance
1096  level of alpha.
1097 
1098  Ref http://www.itl.nist.gov/div898/handbook/ 1.3.6.6.5
1099  param x probability or significance level of the test, >=0 and < 1
1100  param n1 degrees of freedom of first sample, n1 > 0
1101  param n2 degrees of freedom of second sample, n2 > 0
1102  return the statistic (a ratio variance1/variance2) at this prob */
1103  double FDistPDF(double x, int n1, int n2)
1104  {
1105  try
1106  {
1107  double dn1(n1), dn2(n2);
1108  double F =
1109  Gamma((dn1 + dn2) / 2.0) / (Gamma(dn1 / 2.0) * Gamma(dn2 / 2.0));
1110  F *= ::pow(dn1 / dn2, dn1 / 2.0) * ::pow(x, dn1 / 2.0 - 1.0);
1111  F /= ::pow(1.0 + x * dn1 / dn2, (dn1 + dn2) / 2.0);
1112  return F;
1113  }
1114  catch (Exception& e)
1115  {
1116  GNSSTK_RETHROW(e);
1117  }
1118  }
1119 
1120  /* F-distribution percent point function, or inverse of the F-dist CDF.
1121  this function(prob,N1,N2) == F where prob = FDistCDF(F,N1,N2).
1122  Ref http://www.itl.nist.gov/div898/handbook/ 1.3.6.6.5
1123  param prob probability or significance level of the test, >=0 and < 1
1124  param n1 degrees of freedom of first sample, n1 > 0
1125  param n2 degrees of freedom of second sample, n2 > 0
1126  return F the statistic (a ratio variance1/variance2) at this prob */
1127  double invFDistCDF(double prob, int n1, int n2)
1128  {
1129  try
1130  {
1131  if (prob < 0 || prob >= 1)
1132  {
1133  GNSSTK_THROW(Exception("Invalid probability argument"));
1134  }
1135  if (n1 <= 0 || n2 <= 0)
1136  {
1137  GNSSTK_THROW(Exception("Non-positive degree of freedom"));
1138  }
1139 
1140  static const double eps(100000 *
1141  std::numeric_limits<double>().epsilon());
1142  if (prob < eps)
1143  {
1144  return 0.0;
1145  }
1146  if (1.0 - prob < eps)
1147  {
1148  GNSSTK_THROW(Exception("Invalid probability -- too close to 1.0"));
1149  }
1150 
1151  // find F such that prob == FDistCDF(F,n1,n2); use bracket method
1152 
1153  /* NB Abramowitz and Stegan 26.6.9: FDistCDF(F,n1,n2) =
1154  1-FDistCDF(1/F,n2,n1) */
1155  bool swap(false);
1156  int N1(n1), N2(n2);
1157  double alpha(prob);
1158  if (prob < 0.5)
1159  {
1160  swap = true;
1161  N1 = n2, N2 = n1;
1162  alpha = 1.0 - prob;
1163  }
1164 
1165  // we know a0 = FDistCDF(F0,N1,N2) where a0=F0=0 and alpha >= 0.5
1166  double F, F0(0.0), F1, a;
1167  // first find F1 such that a1 = FDistCDF(F1,N1,N2) and a1 > alpha
1168  F1 = 2.0;
1169  while ((a = FDistCDF(F1, N1, N2)) <= alpha)
1170  {
1171  F1 *= 2.0;
1172  }
1173 
1174  // bracket
1175  int n(0); // iteration count
1176  while (1)
1177  {
1178  F = (F0 + F1) / 2.0;
1179  a = FDistCDF(F, N1, N2);
1180  // LOG(INFO) << "LOOP invFDistCDF F = " << n << " " << std::fixed
1181  //<< std::setprecision(15) << F0 << " < " << F << " < " << F1
1182  //<< " and a = " << a << " alpha = " << alpha << " a-alpha = "
1183  //<< std::scientific << std::setprecision(2) << a-alpha << " eps "
1184  //<< eps
1185  //<< " fabs(a-alpha)-eps " << ::fabs(alpha-a)-eps;
1186  if (::fabs(alpha - a) < eps)
1187  {
1188  break;
1189  }
1190  if (a > alpha)
1191  {
1192  F1 = F;
1193  }
1194  else
1195  {
1196  F0 = F;
1197  }
1198  n++;
1199  if (n > 1000)
1200  {
1201  GNSSTK_THROW(Exception("Failed to converge"));
1202  }
1203  }
1204 
1205  return (swap ? 1.0 / F : F);
1206  }
1207  catch (Exception& e)
1208  {
1209  GNSSTK_RETHROW(e);
1210  }
1211  }
1212 
1213 } // namespace gnsstk
gnsstk::invChisqCDF
double invChisqCDF(double alpha, int n)
Definition: SpecialFuncs.cpp:789
gnsstk::errorFunc
double errorFunc(double x)
Definition: SpecialFuncs.cpp:396
gnsstk::invFDistCDF
double invFDistCDF(double prob, int n1, int n2)
Definition: SpecialFuncs.cpp:1127
gnsstk::binomialCoeff
double binomialCoeff(int n, int k)
Definition: SpecialFuncs.cpp:177
gnsstk::ChisqCDF
double ChisqCDF(double x, int n)
Definition: SpecialFuncs.cpp:760
gnsstk::sum
T sum(const ConstVectorBase< T, BaseClass > &l)
Definition: VectorBaseOperators.hpp:84
gnsstk::NormalCDF
double NormalCDF(double x, double mu, double sig)
Definition: SpecialFuncs.cpp:588
gnsstk::lnGamma
double lnGamma(double x)
Definition: SpecialFuncs.cpp:64
gnsstk::incompGamma
double incompGamma(double a, double x)
Definition: SpecialFuncs.cpp:332
gnsstk::cfIBeta
double cfIBeta(double x, double a, double b)
Definition: SpecialFuncs.cpp:434
gnsstk
For Sinex::InputHistory.
Definition: BasicFramework.cpp:50
gnsstk::incompleteBeta
double incompleteBeta(double x, double a, double b)
Definition: SpecialFuncs.cpp:499
gnsstk::Exception
Definition: Exception.hpp:151
gnsstk::factorial
double factorial(int n)
Definition: SpecialFuncs.cpp:114
gnsstk::compIncompGamma
double compIncompGamma(double a, double x)
Definition: SpecialFuncs.cpp:366
gnsstk::seriesIncompGamma
double seriesIncompGamma(double a, double x)
Definition: SpecialFuncs.cpp:222
y
page HOWTO subpage DoxygenGuide Documenting Your Code page DoxygenGuide Documenting Your Code todo Flesh out this document section doctips Tips for Documenting When defining make sure that the prototype is identical between the cpp and hpp including both the namespaces and the parameter names for you have std::string as the return type in the hpp file and string as the return type in the cpp Doxygen may get confused and autolink to the cpp version with no documentation If you don t use the same parameter names between the cpp and hpp that will also confuse Doxygen Don t put type information in return or param documentation It doesn t really add anything and will often cause Doxygen to complain and not produce the documentation< br > use note Do not put a comma after a param name unless you mean to document multiple parameters< br/> the output stream</code >< br/> y
Definition: DOCUMENTING.dox:15
gnsstk::FDistPDF
double FDistPDF(double x, int n1, int n2)
Definition: SpecialFuncs.cpp:1103
gnsstk::FDistCDF
double FDistCDF(double F, int n1, int n2)
Definition: SpecialFuncs.cpp:1050
log
#define log
Definition: DiscCorr.cpp:625
gnsstk::min
T min(const SparseMatrix< T > &SM)
Maximum element - return 0 if empty.
Definition: SparseMatrix.hpp:858
arg
double arg
Definition: IERS1996UT1mUTCData.hpp:48
gnsstk::compErrorFunc
double compErrorFunc(double x)
Definition: SpecialFuncs.cpp:415
gnsstk::lnFactorial
double lnFactorial(int n)
Definition: SpecialFuncs.cpp:149
gnsstk::beta
double beta(double x, double y)
Definition: SpecialFuncs.cpp:204
GNSSTK_RETHROW
#define GNSSTK_RETHROW(exc)
Definition: Exception.hpp:369
gnsstk::invNormalCDF
double invNormalCDF(double prob, double mu, double sig)
Definition: SpecialFuncs.cpp:615
gnsstk::contfracIncompGamma
double contfracIncompGamma(double a, double x)
Definition: SpecialFuncs.cpp:268
gnsstk::ChisqPDF
double ChisqPDF(double x, int n)
Definition: SpecialFuncs.cpp:732
GNSSTK_THROW
#define GNSSTK_THROW(exc)
Definition: Exception.hpp:366
gnsstk::Gamma
double Gamma(double x)
Definition: SpecialFuncs.cpp:98
example3.mu
int mu
Definition: example3.py:36
gnsstk::NormalPDF
double NormalPDF(double x, double mu, double sig)
Definition: SpecialFuncs.cpp:570
gnsstk::StudentsCDF
double StudentsCDF(double t, int n)
Definition: SpecialFuncs.cpp:921
gnsstk::StudentsPDF
double StudentsPDF(double X, int n)
Definition: SpecialFuncs.cpp:890
SpecialFuncs.hpp
gnsstk::invStudentsCDF
double invStudentsCDF(double prob, int n)
Definition: SpecialFuncs.cpp:950

gnsstk
Author(s):
autogenerated on Wed Oct 25 2023 02:40:41