1 /*
2 * Scythe Statistical Library Copyright (C) 2000-2002 Andrew D. Martin
3 * and Kevin M. Quinn; 2002-present Andrew D. Martin, Kevin M. Quinn,
4 * and Daniel Pemstein. All Rights Reserved.
5 *
6 * This program is free software; you can redistribute it and/or
7 * modify under the terms of the GNU General Public License as
8 * published by Free Software Foundation; either version 2 of the
9 * License, or (at your option) any later version. See the text files
10 * COPYING and LICENSE, distributed with this source code, for further
11 * information.
12 * --------------------------------------------------------------------
13 * scythestat/distributions.h
14 *
15 */
16
17 /*! \file distributions.h
18 *
19 * \brief Definitions for probability density functions
20 * (PDFs), cumulative distribution functions (CDFs), and some common
21 * functions (gamma, beta, etc).
22 *
23 * This file provides functions that evaluate the PDFs and CDFs of a
24 * number of probability distributions. In addition, it includes
25 * definitions for another of related functions, such as the gamma
26 * and beta functions.
27 *
28 * The various distribution functions in this file operate on both
29 * scalar quantiles and matrices of quantiles and the
30 * definitions of both forms of these functions appear below. We
31 * provide explicit documentation only for the scalar versions of the
32 * these functions and describe the Matrix versions in the scalar
33 * calls' documents. Much like the operators in matrix.h, we
34 * implement these overloaded versions of the distribution functions
35 * in terms of both generalized and default templates to allow for
36 * explicit control over the template type of the returned Matrix.
37 *
38 * \note Doxygen does not correctly expand the macro definitions we use
39 * to generate the Matrix versions of the various distribution
40 * functions. Therefore, it incorrectly substitutes the macro
41 * variable
42 * __VA_ARGS__ for the actual parameter values in the parameter lists
43 * of each of these functions. For example, the definitions of the
44 * Matrix versions of pbeta are listed as
45 * \code
46 * template<matrix_order RO, matrix_style RS, matrix_order PO, matrix_style PS>
47 * Matrix<double, RO, RS> scythe::pbeta (const Matrix<double, PO, PS> &X, __VA_ARGS__)
48 *
49 * template<matrix_order O, matrix_style S>
50 * Matrix<double, O, Concrete> scythe::pbeta (const Matrix<double, O, S> &X, __VA_ARGS__)
51 * \endcode
52 * when they should be
53 * \code
54 * template<matrix_order RO, matrix_style RS, matrix_order PO, matrix_style PS>
55 * Matrix<double, RO, RS> scythe::pbeta (const Matrix<double, PO, PS> &X, double a, double b)
56 *
57 * template<matrix_order O, matrix_style S>
58 * Matrix<double, O, Concrete> scythe::pbeta (const Matrix<double, O, S> &X, double a, double b)
59 * \endcode
60 *
61 * \par
62 * Furthermore, Doxygen erroneously lists a number of variables at the
63 * end of this document that are not, in fact, declared in
64 * distributions.h. Again, this error is a result of Doxygen's macro
65 * parsing capabilities.
66 *
67 */
68
69 /* TODO: Figure out how to get doxygen to stop listing erroneous
70 * variables at the end of the doc for this file. They stem from it
71 * misreading the nested macro calls used to generate matrix procs.
72 */
73
74 /* TODO: Full R-style versions of these function including arbitrary
75 * recycling of matrix arguments. This is going to have to wait for
76 * variadic templates to be doable without a complete mess. There is
77 * currently a variadic patch available for g++, perhaps we can add a
78 * SCYTHE_VARIADIC flag and include these as option until they become
79 * part of the standard in 2009. Something to come back to after 1.0.
80 */
81
82 #ifndef SCYTHE_DISTRIBUTIONS_H
83 #define SCYTHE_DISTRIBUTIONS_H
84
85 #include <iostream>
86 #include <cmath>
87 #include <cfloat>
88 #include <climits>
89 #include <algorithm>
90 #include <limits>
91
92 #ifdef HAVE_IEEEFP_H
93 #include <ieeefp.h>
94 #endif
95
96 #ifdef SCYTHE_COMPILE_DIRECT
97 #include "matrix.h"
98 #include "ide.h"
99 #include "error.h"
100 #else
101 #include "scythestat/matrix.h"
102 #include "scythestat/ide.h"
103 #include "scythestat/error.h"
104 #endif
105
106 /* Fill in some defs from R that aren't in math.h */
107 #ifndef M_PI
108 #define M_PI 3.141592653589793238462643383280
109 #endif
110 #define M_LN_SQRT_2PI 0.918938533204672741780329736406
111 #define M_LN_SQRT_PId2 0.225791352644727432363097614947
112 #define M_1_SQRT_2PI 0.39894228040143267793994605993
113 #define M_2PI 6.28318530717958647692528676655
114 #define M_SQRT_32 5.656854249492380195206754896838
115
116 #ifndef HAVE_TRUNC
117 /*! @cond */
trunc(double x)118 inline double trunc(double x) throw ()
119 {
120 if (x >= 0)
121 return std::floor(x);
122 else
123 return std::ceil(x);
124 }
125 /*! @endcond */
126 #endif
127
128 /* Many random number generators, pdfs, cdfs, and functions (gamma,
129 * etc) in this file are based on code from the R Project, version
130 * 1.6.0-1.7.1. This code is available under the terms of the GNU
131 * GPL. Original copyright:
132 *
133 * Copyright (C) 1998 Ross Ihaka
134 * Copyright (C) 2000-2002 The R Development Core Team
135 * Copyright (C) 2003 The R Foundation
136 */
137
138 namespace scythe {
139
140 /*! @cond */
141
142 /* Forward declarations */
143 double gammafn (double);
144 double lngammafn (double);
145 double lnbetafn (double, double);
146 double pgamma (double, double, double);
147 double dgamma(double, double, double);
148 double pnorm (double, double, double);
149
150 /*! @endcond */
151
152 /********************
153 * Helper Functions *
154 ********************/
155 namespace {
156
157 /* Evaluate a Chebysheve series at a given point */
158 double
chebyshev_eval(double x,const double * a,int n)159 chebyshev_eval (double x, const double *a, int n)
160 {
161 SCYTHE_CHECK_10(n < 1 || n > 1000, scythe_invalid_arg,
162 "n not on [1, 1000]");
163
164 SCYTHE_CHECK_10(x < -1.1 || x > 1.1, scythe_invalid_arg,
165 "x not on [-1.1, 1.1]");
166
167 double b0, b1, b2;
168 b0 = b1 = b2 = 0;
169
170 double twox = x * 2;
171
172 for (int i = 1; i <= n; ++i) {
173 b2 = b1;
174 b1 = b0;
175 b0 = twox * b1 - b2 + a[n - i];
176 }
177
178 return (b0 - b2) * 0.5;
179 }
180
181 /* Computes the log gamma correction factor for x >= 10 */
182 double
lngammacor(double x)183 lngammacor(double x)
184 {
185 const double algmcs[15] = {
186 +.1666389480451863247205729650822e+0,
187 -.1384948176067563840732986059135e-4,
188 +.9810825646924729426157171547487e-8,
189 -.1809129475572494194263306266719e-10,
190 +.6221098041892605227126015543416e-13,
191 };
192
193 SCYTHE_CHECK_10(x < 10, scythe_invalid_arg,
194 "This function requires x >= 10");
195 SCYTHE_CHECK_10(x >= 3.745194030963158e306, scythe_range_error,
196 "Underflow");
197
198 if (x < 94906265.62425156) {
199 double tmp = 10 / x;
200 return chebyshev_eval(tmp * tmp * 2 - 1, algmcs, 5) / x;
201 }
202
203 return 1 / (x * 12);
204 }
205
206 /* Evaluates the "deviance part" */
207 double
bd0(double x,double np)208 bd0(double x, double np)
209 {
210
211 if(std::fabs(x - np) < 0.1 * (x + np)) {
212 double v = (x - np) / (x + np);
213 double s = (x - np) * v;
214 double ej = 2 * x * v;
215 v = v * v;
216 for (int j = 1; ; j++) {
217 ej *= v;
218 double s1 = s + ej / ((j << 1) + 1);
219 if (s1 == s)
220 return s1;
221 s = s1;
222 }
223 }
224
225 return x * std::log(x / np) + np - x;
226 }
227
228 /* Computes the log of the error term in Stirling's formula */
229 double
stirlerr(double n)230 stirlerr(double n)
231 {
232 #define S0 0.083333333333333333333 /* 1/12 */
233 #define S1 0.00277777777777777777778 /* 1/360 */
234 #define S2 0.00079365079365079365079365 /* 1/1260 */
235 #define S3 0.000595238095238095238095238 /* 1/1680 */
236 #define S4 0.0008417508417508417508417508/* 1/1188 */
237
238 /* error for 0, 0.5, 1.0, 1.5, ..., 14.5, 15.0 */
239 const double sferr_halves[31] = {
240 0.0, /* n=0 - wrong, place holder only */
241 0.1534264097200273452913848, /* 0.5 */
242 0.0810614667953272582196702, /* 1.0 */
243 0.0548141210519176538961390, /* 1.5 */
244 0.0413406959554092940938221, /* 2.0 */
245 0.03316287351993628748511048, /* 2.5 */
246 0.02767792568499833914878929, /* 3.0 */
247 0.02374616365629749597132920, /* 3.5 */
248 0.02079067210376509311152277, /* 4.0 */
249 0.01848845053267318523077934, /* 4.5 */
250 0.01664469118982119216319487, /* 5.0 */
251 0.01513497322191737887351255, /* 5.5 */
252 0.01387612882307074799874573, /* 6.0 */
253 0.01281046524292022692424986, /* 6.5 */
254 0.01189670994589177009505572, /* 7.0 */
255 0.01110455975820691732662991, /* 7.5 */
256 0.010411265261972096497478567, /* 8.0 */
257 0.009799416126158803298389475, /* 8.5 */
258 0.009255462182712732917728637, /* 9.0 */
259 0.008768700134139385462952823, /* 9.5 */
260 0.008330563433362871256469318, /* 10.0 */
261 0.007934114564314020547248100, /* 10.5 */
262 0.007573675487951840794972024, /* 11.0 */
263 0.007244554301320383179543912, /* 11.5 */
264 0.006942840107209529865664152, /* 12.0 */
265 0.006665247032707682442354394, /* 12.5 */
266 0.006408994188004207068439631, /* 13.0 */
267 0.006171712263039457647532867, /* 13.5 */
268 0.005951370112758847735624416, /* 14.0 */
269 0.005746216513010115682023589, /* 14.5 */
270 0.005554733551962801371038690 /* 15.0 */
271 };
272 double nn;
273
274 if (n <= 15.0) {
275 nn = n + n;
276 if (nn == (int)nn)
277 return(sferr_halves[(int)nn]);
278 return (scythe::lngammafn(n + 1.) - (n + 0.5) * std::log(n) + n -
279 std::log(std::sqrt(2 * M_PI)));
280 }
281
282 nn = n*n;
283 if (n > 500)
284 return((S0 - S1 / nn) / n);
285 if (n > 80)
286 return((S0 - (S1 - S2 / nn) / nn) / n);
287 if (n > 35)
288 return((S0 - (S1 - (S2 - S3 / nn) / nn) / nn) / n);
289 /* 15 < n <= 35 : */
290 return((S0 - (S1 - (S2 - (S3 - S4 / nn) / nn) / nn) / nn) / n);
291 #undef S1
292 #undef S2
293 #undef S3
294 #undef S4
295 }
296
297
298 /* Helper for dpois and dgamma */
299 double
dpois_raw(double x,double lambda)300 dpois_raw (double x, double lambda)
301 {
302 if (lambda == 0)
303 return ( (x == 0) ? 1.0 : 0.0);
304
305 if (x == 0)
306 return std::exp(-lambda);
307
308 if (x < 0)
309 return 0.0;
310
311 return std::exp(-stirlerr(x) - bd0(x, lambda))
312 / std::sqrt(2 * M_PI * x);
313 }
314
315
316 /* helper for pbeta */
317 double
pbeta_raw(double x,double pin,double qin)318 pbeta_raw(double x, double pin, double qin)
319 {
320 double ans, c, finsum, p, ps, p1, q, term, xb, xi, y;
321 int n, i, ib, swap_tail;
322
323 const double eps = .5 * DBL_EPSILON;
324 const double sml = DBL_MIN;
325 const double lneps = std::log(eps);
326 const double lnsml = std::log(eps);
327
328 if (pin / (pin + qin) < x) {
329 swap_tail = 1;
330 y = 1 - x;
331 p = qin;
332 q = pin;
333 } else {
334 swap_tail=0;
335 y = x;
336 p = pin;
337 q = qin;
338 }
339
340 if ((p + q) * y / (p + 1) < eps) {
341 ans = 0;
342 xb = p * std::log(std::max(y,sml)) - std::log(p) - lnbetafn(p,q);
343 if (xb > lnsml && y != 0)
344 ans = std::exp(xb);
345 if (swap_tail)
346 ans = 1-ans;
347 } else {
348 ps = q - std::floor(q);
349 if (ps == 0)
350 ps = 1;
351 xb = p * std::log(y) - lnbetafn(ps, p) - std::log(p);
352 ans = 0;
353 if (xb >= lnsml) {
354 ans = std::exp(xb);
355 term = ans * p;
356 if (ps != 1) {
357 n = (int)std::max(lneps/std::log(y), 4.0);
358 for(i = 1; i <= n; i++){
359 xi = i;
360 term *= (xi-ps)*y/xi;
361 ans += term/(p+xi);
362 }
363 }
364 }
365 if (q > 1) {
366 xb = p * std::log(y) + q * std::log(1 - y)
367 - lnbetafn(p, q) - std::log(q);
368 ib = (int) std::max(xb / lnsml, 0.0);
369 term = std::exp(xb - ib * lnsml);
370 c = 1 / (1 - y);
371 p1 = q * c / (p + q - 1);
372
373 finsum = 0;
374 n = (int) q;
375 if(q == n)
376 n--;
377 for (i = 1; i <= n; i++) {
378 if(p1 <= 1 && term / eps <= finsum)
379 break;
380 xi = i;
381 term = (q -xi + 1) * c * term / (p + q - xi);
382 if (term > 1) {
383 ib--;
384 term *= sml;
385 }
386 if (ib == 0)
387 finsum += term;
388 }
389 ans += finsum;
390 }
391
392 if(swap_tail)
393 ans = 1-ans;
394 ans = std::max(std::min(ans,1.),0.);
395 }
396 return ans;
397 }
398
399 /* Helper for dbinom */
400 double
dbinom_raw(double x,double n,double p,double q)401 dbinom_raw (double x, double n, double p, double q)
402 {
403 double f, lc;
404
405 if (p == 0)
406 return((x == 0) ? 1.0 : 0.0);
407 if (q == 0)
408 return((x == n) ? 1.0 : 0.0);
409
410 if (x == 0) {
411 if(n == 0)
412 return 1.0;
413
414 lc = (p < 0.1) ? -bd0(n, n * q) - n * p : n * std::log(q);
415 return(std::exp(lc));
416 }
417 if (x == n) {
418 lc = (q < 0.1) ? -bd0(n,n * p) - n * q : n * std::log(p);
419 return(std::exp(lc));
420 }
421
422 if (x < 0 || x > n)
423 return 0.0;
424
425 lc = stirlerr(n) - stirlerr(x) - stirlerr(n-x) - bd0(x,n*p) -
426 bd0(n - x, n * q);
427
428 f = (M_2PI * x * (n-x)) / n;
429
430 return (std::exp(lc) / std::sqrt(f));
431 }
432
433 /* The normal probability density function implementation. */
434
435 #define SIXTEN 16
436 #define do_del(X) \
437 xsq = trunc(X * SIXTEN) / SIXTEN; \
438 del = (X - xsq) * (X + xsq); \
439 if(log_p) { \
440 *cum = (-xsq * xsq * 0.5) + (-del * 0.5) + std::log(temp); \
441 if((lower && x > 0.) || (upper && x <= 0.)) \
442 *ccum = ::log1p(-std::exp(-xsq * xsq * 0.5) * \
443 std::exp(-del * 0.5) * temp); \
444 } \
445 else { \
446 *cum = std::exp(-xsq * xsq * 0.5) * std::exp(-del * 0.5) * temp; \
447 *ccum = 1.0 - *cum; \
448 }
449
450 #define swap_tail \
451 if (x > 0.) {/* swap ccum <--> cum */ \
452 temp = *cum; if(lower) *cum = *ccum; *ccum = temp; \
453 }
454
455 void
pnorm_both(double x,double * cum,double * ccum,int i_tail,bool log_p)456 pnorm_both(double x, double *cum, double *ccum, int i_tail,
457 bool log_p)
458 {
459 const double a[5] = {
460 2.2352520354606839287,
461 161.02823106855587881,
462 1067.6894854603709582,
463 18154.981253343561249,
464 0.065682337918207449113
465 };
466 const double b[4] = {
467 47.20258190468824187,
468 976.09855173777669322,
469 10260.932208618978205,
470 45507.789335026729956
471 };
472 const double c[9] = {
473 0.39894151208813466764,
474 8.8831497943883759412,
475 93.506656132177855979,
476 597.27027639480026226,
477 2494.5375852903726711,
478 6848.1904505362823326,
479 11602.651437647350124,
480 9842.7148383839780218,
481 1.0765576773720192317e-8
482 };
483 const double d[8] = {
484 22.266688044328115691,
485 235.38790178262499861,
486 1519.377599407554805,
487 6485.558298266760755,
488 18615.571640885098091,
489 34900.952721145977266,
490 38912.003286093271411,
491 19685.429676859990727
492 };
493 const double p[6] = {
494 0.21589853405795699,
495 0.1274011611602473639,
496 0.022235277870649807,
497 0.001421619193227893466,
498 2.9112874951168792e-5,
499 0.02307344176494017303
500 };
501 const double q[5] = {
502 1.28426009614491121,
503 0.468238212480865118,
504 0.0659881378689285515,
505 0.00378239633202758244,
506 7.29751555083966205e-5
507 };
508
509 double xden, xnum, temp, del, eps, xsq, y;
510 int i, lower, upper;
511
512 /* Consider changing these : */
513 eps = DBL_EPSILON * 0.5;
514
515 /* i_tail in {0,1,2} =^= {lower, upper, both} */
516 lower = i_tail != 1;
517 upper = i_tail != 0;
518
519 y = std::fabs(x);
520 if (y <= 0.67448975) {
521 /* qnorm(3/4) = .6744.... -- earlier had 0.66291 */
522 if (y > eps) {
523 xsq = x * x;
524 xnum = a[4] * xsq;
525 xden = xsq;
526 for (i = 0; i < 3; ++i) {
527 xnum = (xnum + a[i]) * xsq;
528 xden = (xden + b[i]) * xsq;
529 }
530 } else xnum = xden = 0.0;
531
532 temp = x * (xnum + a[3]) / (xden + b[3]);
533 if(lower) *cum = 0.5 + temp;
534 if(upper) *ccum = 0.5 - temp;
535 if(log_p) {
536 if(lower) *cum = std::log(*cum);
537 if(upper) *ccum = std::log(*ccum);
538 }
539 } else if (y <= M_SQRT_32) {
540 /* Evaluate pnorm for 0.674.. = qnorm(3/4) < |x| <= sqrt(32)
541 * ~= 5.657 */
542
543 xnum = c[8] * y;
544 xden = y;
545 for (i = 0; i < 7; ++i) {
546 xnum = (xnum + c[i]) * y;
547 xden = (xden + d[i]) * y;
548 }
549 temp = (xnum + c[7]) / (xden + d[7]);
550 do_del(y);
551 swap_tail;
552 } else if (log_p
553 || (lower && -37.5193 < x && x < 8.2924)
554 || (upper && -8.2929 < x && x < 37.5193)
555 ) {
556 /* Evaluate pnorm for x in (-37.5, -5.657) union (5.657, 37.5) */
557 xsq = 1.0 / (x * x);
558 xnum = p[5] * xsq;
559 xden = xsq;
560 for (i = 0; i < 4; ++i) {
561 xnum = (xnum + p[i]) * xsq;
562 xden = (xden + q[i]) * xsq;
563 }
564 temp = xsq * (xnum + p[4]) / (xden + q[4]);
565 temp = (M_1_SQRT_2PI - temp) / y;
566 do_del(x);
567 swap_tail;
568 } else {
569 if (x > 0) {
570 *cum = 1.;
571 *ccum = 0.;
572 } else {
573 *cum = 0.;
574 *ccum = 1.;
575 }
576 //XXX commented out for debug-on testing of ordfactanal
577 //(and perhaps others) since they tend to throw on the first
578 //iteration
579 //SCYTHE_THROW_10(scythe_convergence_error, "Did not converge");
580 }
581
582 return;
583 }
584 #undef SIXTEN
585 #undef do_del
586 #undef swap_tail
587
588 /* The standard normal distribution function */
589 double
pnorm1(double x,bool lower_tail,bool log_p)590 pnorm1 (double x, bool lower_tail, bool log_p)
591 {
592 SCYTHE_CHECK_10(! R_finite(x), scythe_invalid_arg,
593 "Quantile x is inifinte (+/-Inf) or NaN");
594
595 double p, cp;
596 pnorm_both(x, &p, &cp, (lower_tail ? 0 : 1), log_p);
597
598 return (lower_tail ? p : cp);
599 }
600 } // anonymous namespace
601
602 /*************
603 * Functions *
604 *************/
605
606 /* The gamma function */
607
608 /*! \brief The gamma function.
609 *
610 * Computes the gamma function, evaluated at \a x.
611 *
612 * \param x The value to compute gamma at.
613 *
614 * \see lngammafn(double x)
615 * \see pgamma(double x, double shape, double scale)
616 * \see dgamma(double x, double shape, double scale)
617 * \see rng::rgamma(double shape, double scale)
618 *
619 * \throw scythe_range_error (Level 1)
620 * \throw scythe_precision_error (Level 1)
621 */
622 inline double
gammafn(double x)623 gammafn (double x)
624 {
625 const double gamcs[22] = {
626 +.8571195590989331421920062399942e-2,
627 +.4415381324841006757191315771652e-2,
628 +.5685043681599363378632664588789e-1,
629 -.4219835396418560501012500186624e-2,
630 +.1326808181212460220584006796352e-2,
631 -.1893024529798880432523947023886e-3,
632 +.3606925327441245256578082217225e-4,
633 -.6056761904460864218485548290365e-5,
634 +.1055829546302283344731823509093e-5,
635 -.1811967365542384048291855891166e-6,
636 +.3117724964715322277790254593169e-7,
637 -.5354219639019687140874081024347e-8,
638 +.9193275519859588946887786825940e-9,
639 -.1577941280288339761767423273953e-9,
640 +.2707980622934954543266540433089e-10,
641 -.4646818653825730144081661058933e-11,
642 +.7973350192007419656460767175359e-12,
643 -.1368078209830916025799499172309e-12,
644 +.2347319486563800657233471771688e-13,
645 -.4027432614949066932766570534699e-14,
646 +.6910051747372100912138336975257e-15,
647 -.1185584500221992907052387126192e-15,
648 };
649
650
651 double y = std::fabs(x);
652
653 if (y <= 10) {
654
655 /* Compute gamma(x) for -10 <= x <= 10
656 * Reduce the interval and find gamma(1 + y) for 0 <= y < 1
657 * first of all. */
658
659 int n = (int) x;
660 if (x < 0)
661 --n;
662
663 y = x - n;/* n = floor(x) ==> y in [ 0, 1 ) */
664 --n;
665 double value = chebyshev_eval(y * 2 - 1, gamcs, 22) + .9375;
666
667 if (n == 0)
668 return value;/* x = 1.dddd = 1+y */
669
670 if (n < 0) {
671 /* compute gamma(x) for -10 <= x < 1 */
672
673 /* If the argument is exactly zero or a negative integer */
674 /* then return NaN. */
675 SCYTHE_CHECK_10(x == 0 || (x < 0 && x == n + 2),
676 scythe_range_error, "x is 0 or a negative integer");
677
678 /* The answer is less than half precision */
679 /* because x too near a negative integer. */
680 SCYTHE_CHECK_10(x < -0.5 &&
681 std::fabs(x - (int)(x - 0.5) / x) < 67108864.0,
682 scythe_precision_error,
683 "Answer < 1/2 precision because x is too near" <<
684 " a negative integer");
685
686 /* The argument is so close to 0 that the result
687 * * would overflow. */
688 SCYTHE_CHECK_10(y < 2.2474362225598545e-308, scythe_range_error,
689 "x too close to 0");
690
691 n = -n;
692
693 for (int i = 0; i < n; i++)
694 value /= (x + i);
695
696 return value;
697 } else {
698 /* gamma(x) for 2 <= x <= 10 */
699
700 for (int i = 1; i <= n; i++) {
701 value *= (y + i);
702 }
703 return value;
704 }
705 } else {
706 /* gamma(x) for y = |x| > 10. */
707
708 /* Overflow */
709 SCYTHE_CHECK_10(x > 171.61447887182298,
710 scythe_range_error,"Overflow");
711
712 /* Underflow */
713 SCYTHE_CHECK_10(x < -170.5674972726612,
714 scythe_range_error, "Underflow");
715
716 double value = std::exp((y - 0.5) * std::log(y) - y
717 + M_LN_SQRT_2PI + lngammacor(y));
718
719 if (x > 0)
720 return value;
721
722 SCYTHE_CHECK_10(std::fabs((x - (int)(x - 0.5))/x) < 67108864.0,
723 scythe_precision_error,
724 "Answer < 1/2 precision because x is " <<
725 "too near a negative integer");
726
727 double sinpiy = std::sin(M_PI * y);
728
729 /* Negative integer arg - overflow */
730 SCYTHE_CHECK_10(sinpiy == 0, scythe_range_error, "Overflow");
731
732 return -M_PI / (y * sinpiy * value);
733 }
734 }
735
736 /* The natural log of the absolute value of the gamma function */
737 /*! \brief The natural log of the absolute value of the gamma
738 * function.
739 *
740 * Computes the natural log of the absolute value of the gamma
741 * function, evaluated at \a x.
742 *
743 * \param x The value to compute log(abs(gamma())) at.
744 *
745 * \see gammafn(double x)
746 * \see pgamma(double x, double shape, double scale)
747 * \see dgamma(double x, double shape, double scale)
748 * \see rng::rgamma(double shape, double scale)
749 *
750 * \throw scythe_range_error (Level 1)
751 * \throw scythe_precision_error (Level 1)
752 */
753 inline double
lngammafn(double x)754 lngammafn(double x)
755 {
756 SCYTHE_CHECK_10(x <= 0 && x == (int) x, scythe_range_error,
757 "x is 0 or a negative integer");
758
759 double y = std::fabs(x);
760
761 if (y <= 10)
762 return std::log(std::fabs(gammafn(x)));
763
764 SCYTHE_CHECK_10(y > 2.5327372760800758e+305, scythe_range_error,
765 "Overflow");
766
767 if (x > 0) /* i.e. y = x > 10 */
768 return M_LN_SQRT_2PI + (x - 0.5) * std::log(x) - x
769 + lngammacor(x);
770
771 /* else: x < -10; y = -x */
772 double sinpiy = std::fabs(std::sin(M_PI * y));
773
774 if (sinpiy == 0) /* Negative integer argument */
775 throw scythe_exception("UNEXPECTED ERROR",
776 __FILE__, __func__, __LINE__,
777 "ERROR: Should never happen!");
778
779 double ans = M_LN_SQRT_PId2 + (x - 0.5) * std::log(y) - x - std::log(sinpiy)
780 - lngammacor(y);
781
782 SCYTHE_CHECK_10(std::fabs((x - (int)(x - 0.5)) * ans / x)
783 < 1.490116119384765696e-8, scythe_precision_error,
784 "Answer < 1/2 precision because x is "
785 << "too near a negative integer");
786
787 return ans;
788 }
789
790 /* The beta function */
791 /*! \brief The beta function.
792 *
793 * Computes beta function, evaluated at (\a a, \a b).
794 *
795 * \param a The first parameter.
796 * \param b The second parameter.
797 *
798 * \see lnbetafn(double a, double b)
799 * \see pbeta(double x, double a, double b)
800 * \see dbeta(double x, double a, double b)
801 * \see rng::rbeta(double a, double b)
802 *
803 * \throw scythe_invalid_arg (Level 1)
804 * \throw scythe_range_error (Level 1)
805 * \throw scythe_precision_error (Level 1)
806 */
807 inline double
betafn(double a,double b)808 betafn(double a, double b)
809 {
810 SCYTHE_CHECK_10(a <= 0 || b <= 0, scythe_invalid_arg, "a or b < 0");
811
812 if (a + b < 171.61447887182298) /* ~= 171.61 for IEEE */
813 return gammafn(a) * gammafn(b) / gammafn(a+b);
814
815 double val = lnbetafn(a, b);
816 SCYTHE_CHECK_10(val < -708.39641853226412, scythe_range_error,
817 "Underflow");
818
819 return std::exp(val);
820 }
821
822 /* The natural log of the beta function */
823 /*! \brief The natural log of the beta function.
824 *
825 * Computes the natural log of the beta function,
826 * evaluated at (\a a, \a b).
827 *
828 * \param a The first parameter.
829 * \param b The second parameter.
830 *
831 * \see betafn(double a, double b)
832 * \see pbeta(double x, double a, double b)
833 * \see dbeta(double x, double a, double b)
834 * \see rng::rbeta(double a, double b)
835 *
836 * \throw scythe_invalid_arg (Level 1)
837 * \throw scythe_range_error (Level 1)
838 * \throw scythe_precision_error (Level 1)
839 */
840 inline double
lnbetafn(double a,double b)841 lnbetafn (double a, double b)
842 {
843 double p, q;
844
845 p = q = a;
846 if(b < p) p = b;/* := min(a,b) */
847 if(b > q) q = b;/* := max(a,b) */
848
849 SCYTHE_CHECK_10(p <= 0 || q <= 0,scythe_invalid_arg, "a or b <= 0");
850
851 if (p >= 10) {
852 /* p and q are big. */
853 double corr = lngammacor(p) + lngammacor(q) - lngammacor(p + q);
854 return std::log(q) * -0.5 + M_LN_SQRT_2PI + corr
855 + (p - 0.5) * std::log(p / (p + q)) + q * std::log(1 + (-p / (p + q)));
856 } else if (q >= 10) {
857 /* p is small, but q is big. */
858 double corr = lngammacor(q) - lngammacor(p + q);
859 return lngammafn(p) + corr + p - p * std::log(p + q)
860 + (q - 0.5) * std::log(1 + (-p / (p + q)));
861 }
862
863 /* p and q are small: p <= q > 10. */
864 return std::log(gammafn(p) * (gammafn(q) / gammafn(p + q)));
865 }
866
867 /* Compute the factorial of a non-negative integer */
868 /*! \brief The factorial function.
869 *
870 * Computes the factorial of \a n.
871 *
872 * \param n The non-negative integer value to compute the factorial of.
873 *
874 * \see lnfactorial(unsigned int n)
875 *
876 */
877 inline int
factorial(unsigned int n)878 factorial (unsigned int n)
879 {
880 if (n == 0)
881 return 1;
882
883 return n * factorial(n - 1);
884 }
885
886 /* Compute the natural log of the factorial of a non-negative
887 * integer
888 */
889 /*! \brief The log of the factorial function.
890 *
891 * Computes the natural log of the factorial of \a n.
892 *
893 * \param n The non-negative integer value to compute the natural log of the factorial of.
894 *
895 * \see factorial(unsigned int n)
896 *
897 */
898 inline double
lnfactorial(unsigned int n)899 lnfactorial (unsigned int n)
900 {
901 double x = n+1;
902 double cof[6] = {
903 76.18009172947146, -86.50532032941677,
904 24.01409824083091, -1.231739572450155,
905 0.1208650973866179e-2, -0.5395239384953e-5
906 };
907 double y = x;
908 double tmp = x + 5.5 - (x + 0.5) * std::log(x + 5.5);
909 double ser = 1.000000000190015;
910 for (int j = 0; j <= 5; j++) {
911 ser += (cof[j] / ++y);
912 }
913 return(std::log(2.5066282746310005 * ser / x) - tmp);
914 }
915
916 /*********************************
917 * Fully Specified Distributions *
918 *********************************/
919
920 /* These macros provide a nice shorthand for the matrix versions of
921 * the pdf and cdf functions.
922 */
923
924 #define SCYTHE_ARGSET(...) __VA_ARGS__
925
926 #define SCYTHE_DISTFUN_MATRIX(NAME, XTYPE, ARGNAMES, ...) \
927 template <matrix_order RO, matrix_style RS, \
928 matrix_order PO, matrix_style PS> \
929 Matrix<double, RO, RS> \
930 NAME (const Matrix<XTYPE, PO, PS>& X, __VA_ARGS__) \
931 { \
932 Matrix<double, RO, Concrete> ret(X.rows(), X.cols(), false); \
933 const_matrix_forward_iterator<XTYPE,RO,PO,PS> xit; \
934 const_matrix_forward_iterator<XTYPE,RO,PO,PS> xlast \
935 = X.template end_f<RO>(); \
936 typename Matrix<double,RO,Concrete>::forward_iterator rit \
937 = ret.begin_f(); \
938 for (xit = X.template begin_f<RO>(); xit != xlast; ++xit) { \
939 *rit = NAME (*xit, ARGNAMES); \
940 ++rit; \
941 } \
942 SCYTHE_VIEW_RETURN(double, RO, RS, ret) \
943 } \
944 \
945 template <matrix_order O, matrix_style S> \
946 Matrix<double, O, Concrete> \
947 NAME (const Matrix<XTYPE, O, S>& X, __VA_ARGS__) \
948 { \
949 return NAME <O, Concrete, O, S> (X, ARGNAMES); \
950 }
951
952 /**** The Beta Distribution ****/
953
954 /* CDFs */
955
956 /*! \brief The beta distribution function.
957 *
958 * Computes the value of the beta cumulative distribution function
959 * with shape parameters \a a and \a b at the desired quantile,
960 * \a x.
961 *
962 * It is also possible to call this function with a Matrix of
963 * doubles as its first argument. In this case the function will
964 * return a Matrix of doubles of the same dimension as \a x,
965 * containing the result of evaluating this function at each value
966 * in \a x, given the remaining fixed parameters. By default, the
967 * returned Matrix will be concrete and have the same matrix_order
968 * as \a x, but you may invoke a generalized version of the function
969 * with an explicit template call.
970 *
971 * \param x The desired quantile, between 0 and 1.
972 * \param a The first non-negative beta shape parameter.
973 * \param b The second non-negative beta shape parameter.
974 *
975 * \see dbeta(double x, double a, double b)
976 * \see rng::rbeta(double a, double b)
977 * \see betafn(double a, double b)
978 * \see lnbetafn(double a, double b)
979 *
980 * \throw scythe_invalid_arg (Level 1)
981 * \throw scythe_range_error (Level 1)
982 * \throw scythe_precision_error (Level 1)
983 */
984 inline double
pbeta(double x,double a,double b)985 pbeta(double x, double a, double b)
986 {
987 SCYTHE_CHECK_10(a <= 0 || b <= 0,scythe_invalid_arg, "a or b <= 0");
988
989 if (x <= 0)
990 return 0.;
991 if (x >= 1)
992 return 1.;
993
994 return pbeta_raw(x,a,b);
995 }
996
SCYTHE_DISTFUN_MATRIX(pbeta,double,SCYTHE_ARGSET (a,b),double a,double b)997 SCYTHE_DISTFUN_MATRIX(pbeta, double, SCYTHE_ARGSET(a, b), double a, double b)
998
999 /* PDFs */
1000 /*! \brief The beta density function.
1001 *
1002 * Computes the value of the beta probability density function
1003 * with shape parameters \a a and \a b at the desired quantile,
1004 * \a x.
1005 *
1006 * It is also possible to call this function with a Matrix of
1007 * doubles as its first argument. In this case the function will
1008 * return a Matrix of doubles of the same dimension as \a x,
1009 * containing the result of evaluating this function at each value
1010 * in \a x, given the remaining fixed parameters. By default, the
1011 * returned Matrix will be concrete and have the same matrix_order
1012 * as \a x, but you may invoke a generalized version of the function
1013 * with an explicit template call.
1014 *
1015 * \param x The desired quantile, between 0 and 1.
1016 * \param a The first non-negative beta shape parameter.
1017 * \param b The second non-negative beta shape parameter.
1018 *
1019 * \see pbeta(double x, double a, double b)
1020 * \see rng::rbeta(double a, double b)
1021 * \see betafn(double a, double b)
1022 * \see lnbetafn(double a, double b)
1023 *
1024 * \throw scythe_invalid_arg (Level 1)
1025 * \throw scythe_range_error (Level 1)
1026 * \throw scythe_precision_error (Level 1)
1027 */
1028 inline double
1029 dbeta(double x, double a, double b)
1030 {
1031 SCYTHE_CHECK_10((x < 0.0) || (x > 1.0), scythe_invalid_arg,
1032 "x not in [0,1]");
1033 SCYTHE_CHECK_10(a < 0.0, scythe_invalid_arg, "a < 0");
1034 SCYTHE_CHECK_10(b < 0.0, scythe_invalid_arg, "b < 0");
1035
1036 return (std::pow(x, (a-1.0)) * std::pow((1.0-x), (b-1.0)) )
1037 / betafn(a,b);
1038 }
1039
SCYTHE_DISTFUN_MATRIX(dbeta,double,SCYTHE_ARGSET (a,b),double a,double b)1040 SCYTHE_DISTFUN_MATRIX(dbeta, double, SCYTHE_ARGSET(a, b), double a, double b)
1041
1042 /* Returns the natural log of the ordinate of the Beta density
1043 * evaluated at x with Shape1 a, and Shape2 b
1044 */
1045
1046 /*! \brief The natural log of the ordinate of the beta density
1047 * function.
1048 *
1049 * Computes the value of the natural log of the ordinate of the beta
1050 * probability density function
1051 * with shape parameters \a a and \a b at the desired quantile,
1052 * \a x.
1053 *
1054 * It is also possible to call this function with a Matrix of
1055 * doubles as its first argument. In this case the function will
1056 * return a Matrix of doubles of the same dimension as \a x,
1057 * containing the result of evaluating this function at each value
1058 * in \a x, given the remaining fixed parameters. By default, the
1059 * returned Matrix will be concrete and have the same matrix_order
1060 * as \a x, but you may invoke a generalized version of the function
1061 * with an explicit template call.
1062 *
1063 * \param x The desired quantile, between 0 and 1.
1064 * \param a The first non-negative beta shape parameter.
1065 * \param b The second non-negative beta shape parameter.
1066 *
1067 * \see dbeta(double x, double a, double b)
1068 *
1069 * \throw scythe_invalid_arg (Level 1)
1070 * \throw scythe_range_error (Level 1)
1071 * \throw scythe_precision_error (Level 1)
1072 */
1073 inline double
1074 lndbeta1(double x, double a, double b)
1075 {
1076 SCYTHE_CHECK_10((x < 0.0) || (x > 1.0), scythe_invalid_arg,
1077 "x not in [0,1]");
1078 SCYTHE_CHECK_10(a < 0.0, scythe_invalid_arg, "a < 0");
1079 SCYTHE_CHECK_10(b < 0.0, scythe_invalid_arg, "b < 0");
1080
1081 return (a-1.0) * std::log(x) + (b-1) * std::log(1.0-x)
1082 - lnbetafn(a,b);
1083 }
1084
SCYTHE_DISTFUN_MATRIX(lndbeta1,double,SCYTHE_ARGSET (a,b),double a,double b)1085 SCYTHE_DISTFUN_MATRIX(lndbeta1, double, SCYTHE_ARGSET(a, b), double a, double b)
1086
1087
1088 /**** The Binomial Distribution ****/
1089
1090 /* CDFs */
1091
1092 /*! \brief The binomial distribution function.
1093 *
1094 * Computes the value of the binomial cumulative distribution function
1095 * with \a n trials and \a p probability of success on each trial,
1096 * at the desired quantile \a x.
1097 *
1098 * It is also possible to call this function with a Matrix of
1099 * doubles as its first argument. In this case the function will
1100 * return a Matrix of doubles of the same dimension as \a x,
1101 * containing the result of evaluating this function at each value
1102 * in \a x, given the remaining fixed parameters. By default, the
1103 * returned Matrix will be concrete and have the same matrix_order
1104 * as \a x, but you may invoke a generalized version of the function
1105 * with an explicit template call.
1106 *
1107 * \param x The desired quantile.
1108 * \param n The number of trials.
1109 * \param p The probability of success on each trial.
1110 *
1111 * \see dbinom(double x, unsigned int n, double p)
1112 * \see rng::rbinom(unsigned int n, double p)
1113 *
1114 * \throw scythe_invalid_arg (Level 1)
1115 * \throw scythe_range_error (Level 1)
1116 * \throw scythe_precision_error (Level 1)
1117 */
1118 inline double
1119 pbinom(double x, unsigned int n, double p)
1120 {
1121
1122 SCYTHE_CHECK_10(p < 0 || p > 1, scythe_invalid_arg, "p not in [0,1]");
1123 double X = std::floor(x);
1124
1125 if (X < 0.0)
1126 return 0;
1127
1128 if (n <= X)
1129 return 1;
1130
1131 return pbeta(1 - p, n - X, X + 1);
1132 }
1133
SCYTHE_DISTFUN_MATRIX(pbinom,double,SCYTHE_ARGSET (n,p),unsigned int n,double p)1134 SCYTHE_DISTFUN_MATRIX(pbinom, double, SCYTHE_ARGSET(n, p), unsigned int n, double p)
1135
1136 /* PDFs */
1137 /*! \brief The binomial density function.
1138 *
1139 * Computes the value of the binomial probability density function
1140 * with \a n trials and \a p probability of success on each trial,
1141 * at the desired quantile \a x.
1142 *
1143 * It is also possible to call this function with a Matrix of
1144 * doubles as its first argument. In this case the function will
1145 * return a Matrix of doubles of the same dimension as \a x,
1146 * containing the result of evaluating this function at each value
1147 * in \a x, given the remaining fixed parameters. By default, the
1148 * returned Matrix will be concrete and have the same matrix_order
1149 * as \a x, but you may invoke a generalized version of the function
1150 * with an explicit template call.
1151 *
1152 * \param x The desired quantile.
1153 * \param n The number of trials.
1154 * \param p The probability of success on each trial.
1155 *
1156 * \see pbinom(double x, unsigned int n, double p)
1157 * \see rng::rbinom(unsigned int n, double p)
1158 *
1159 * \throw scythe_invalid_arg (Level 1)
1160 * \throw scythe_range_error (Level 1)
1161 * \throw scythe_precision_error (Level 1)
1162 */
1163 inline double
1164 dbinom(double x, unsigned int n, double p)
1165 {
1166 SCYTHE_CHECK_10(p < 0 || p > 1, scythe_invalid_arg, "p not in [0, 1]");
1167 double X = std::floor(x + 0.5);
1168 return dbinom_raw(X, n, p, 1 - p);
1169 }
1170
SCYTHE_DISTFUN_MATRIX(dbinom,double,SCYTHE_ARGSET (n,p),unsigned int n,double p)1171 SCYTHE_DISTFUN_MATRIX(dbinom, double, SCYTHE_ARGSET(n, p), unsigned int n, double p)
1172
1173 /**** The Chi Squared Distribution ****/
1174
1175 /* CDFs */
1176 /*! \brief The \f$\chi^2\f$ distribution function.
1177 *
1178 * Computes the value of the \f$\chi^2\f$ cumulative distribution
1179 * function with \a df degrees of freedom, at the desired quantile
1180 * \a x.
1181 *
1182 * It is also possible to call this function with a Matrix of
1183 * doubles as its first argument. In this case the function will
1184 * return a Matrix of doubles of the same dimension as \a x,
1185 * containing the result of evaluating this function at each value
1186 * in \a x, given the remaining fixed parameters. By default, the
1187 * returned Matrix will be concrete and have the same matrix_order
1188 * as \a x, but you may invoke a generalized version of the function
1189 * with an explicit template call.
1190 *
1191 * \param x The desired quantile.
1192 * \param df The degrees of freedom.
1193
1194 * \see dchisq(double x, double df)
1195 * \see rng::rchisq(double df)
1196 *
1197 * \throw scythe_invalid_arg (Level 1)
1198 * \throw scythe_range_error (Level 1)
1199 * \throw scythe_precision_error (Level 1)
1200 * \throw scythe_convergence_error (Level 1)
1201 *
1202 */
1203 inline double
1204 pchisq(double x, double df)
1205 {
1206 return pgamma(x, df/2.0, 2.0);
1207 }
1208
SCYTHE_DISTFUN_MATRIX(pchisq,double,df,double df)1209 SCYTHE_DISTFUN_MATRIX(pchisq, double, df, double df)
1210
1211 /* PDFs */
1212 /*! \brief The \f$\chi^2\f$ density function.
1213 *
1214 * Computes the value of the \f$\chi^2\f$ probability density
1215 * function with \a df degrees of freedom, at the desired quantile
1216 * \a x.
1217 *
1218 * It is also possible to call this function with a Matrix of
1219 * doubles as its first argument. In this case the function will
1220 * return a Matrix of doubles of the same dimension as \a x,
1221 * containing the result of evaluating this function at each value
1222 * in \a x, given the remaining fixed parameters. By default, the
1223 * returned Matrix will be concrete and have the same matrix_order
1224 * as \a x, but you may invoke a generalized version of the function
1225 * with an explicit template call.
1226 *
1227 * \param x The desired quantile.
1228 * \param df The degrees of freedom.
1229
1230 * \see pchisq(double x, double df)
1231 * \see rng::rchisq(double df)
1232 *
1233 * \throw scythe_invalid_arg (Level 1)
1234 * \throw scythe_range_error (Level 1)
1235 * \throw scythe_precision_error (Level 1)
1236 * \throw scythe_convergence_error (Level 1)
1237 *
1238 */
1239 inline double
1240 dchisq(double x, double df)
1241 {
1242 return dgamma(x, df / 2.0, 2.0);
1243 }
1244
SCYTHE_DISTFUN_MATRIX(dchisq,double,df,double df)1245 SCYTHE_DISTFUN_MATRIX(dchisq, double, df, double df)
1246
1247 /**** The Exponential Distribution ****/
1248
1249 /* CDFs */
1250 /*! \brief The exponential distribution function.
1251 *
1252 * Computes the value of the exponential cumulative distribution
1253 * function with given \a scale, at the desired quantile
1254 * \a x.
1255 *
1256 * It is also possible to call this function with a Matrix of
1257 * doubles as its first argument. In this case the function will
1258 * return a Matrix of doubles of the same dimension as \a x,
1259 * containing the result of evaluating this function at each value
1260 * in \a x, given the remaining fixed parameters. By default, the
1261 * returned Matrix will be concrete and have the same matrix_order
1262 * as \a x, but you may invoke a generalized version of the function
1263 * with an explicit template call.
1264 *
1265 * \param x The desired quantile.
1266 * \param scale The positive scale of the function.
1267 *
1268 * \see dexp(double x, double scale)
1269 * \see rng::rexp(double scale)
1270 *
1271 * \throw scythe_invalid_arg (Level 1)
1272 */
1273 inline double
1274 pexp(double x, double scale)
1275 {
1276 SCYTHE_CHECK_10(scale <= 0, scythe_invalid_arg, "scale <= 0");
1277
1278 if (x <= 0)
1279 return 0;
1280
1281 return (1 - std::exp(-x*scale));
1282 }
1283
SCYTHE_DISTFUN_MATRIX(pexp,double,scale,double scale)1284 SCYTHE_DISTFUN_MATRIX(pexp, double, scale, double scale)
1285
1286 /* PDFs */
1287 /*! \brief The exponential density function.
1288 *
1289 * Computes the value of the exponential probability density
1290 * function with given \a scale, at the desired quantile
1291 * \a x.
1292 *
1293 * It is also possible to call this function with a Matrix of
1294 * doubles as its first argument. In this case the function will
1295 * return a Matrix of doubles of the same dimension as \a x,
1296 * containing the result of evaluating this function at each value
1297 * in \a x, given the remaining fixed parameters. By default, the
1298 * returned Matrix will be concrete and have the same matrix_order
1299 * as \a x, but you may invoke a generalized version of the function
1300 * with an explicit template call.
1301 *
1302 * \param x The desired quantile.
1303 * \param scale The positive scale of the function.
1304 *
1305 * \see pexp(double x, double scale)
1306 * \see rng::rexp(double scale)
1307 *
1308 * \throw scythe_invalid_arg (Level 1)
1309 */
1310 inline double
1311 dexp(double x, double scale)
1312 {
1313 SCYTHE_CHECK_10(scale <= 0, scythe_invalid_arg, "scale <= 0");
1314
1315 if (x < 0)
1316 return 0;
1317
1318 return std::exp(-x * scale) * scale;
1319 }
1320
SCYTHE_DISTFUN_MATRIX(dexp,double,scale,double scale)1321 SCYTHE_DISTFUN_MATRIX(dexp, double, scale, double scale)
1322
1323 /**** The f Distribution ****/
1324
1325 /* CDFs */
1326 /*! \brief The F distribution function.
1327 *
1328 * Computes the value of the F cumulative distribution function with
1329 * \a df1 and \a df2 degrees of freedom, at the desired quantile \a
1330 * x.
1331 *
1332 * It is also possible to call this function with a Matrix of
1333 * doubles as its first argument. In this case the function will
1334 * return a Matrix of doubles of the same dimension as \a x,
1335 * containing the result of evaluating this function at each value
1336 * in \a x, given the remaining fixed parameters. By default, the
1337 * returned Matrix will be concrete and have the same matrix_order
1338 * as \a x, but you may invoke a generalized version of the function
1339 * with an explicit template call.
1340 *
1341 * \param x The desired quantile.
1342 * \param df1 The non-negative degrees of freedom for the
1343 * \f$\chi^2\f$ variate in the nominator of the F statistic.
1344 * \param df2 The non-negative degrees of freedom for the
1345 * \f$\chi^2\f$ variate in the denominator of the F statistic.
1346 *
1347 *
1348 * \see df(double x, double df1, double df2)
1349 * \see rng::rf(double df1, double df2)
1350 *
1351 * \throw scythe_invalid_arg (Level 1)
1352 * \throw scythe_range_error (Level 1)
1353 * \throw scythe_precision_error (Level 1)
1354 * \throw scythe_convergence_error (Level 1)
1355 */
1356 inline double
1357 pf(double x, double df1, double df2)
1358 {
1359 SCYTHE_CHECK_10(df1 <= 0 || df2 <= 0, scythe_invalid_arg,
1360 "df1 or df2 <= 0");
1361
1362 if (x <= 0)
1363 return 0;
1364
1365 if (df2 > 4e5)
1366 return pchisq(x*df1,df1);
1367 if (df1 > 4e5)
1368 return 1-pchisq(df2/x,df2);
1369
1370 return (1-pbeta(df2 / (df2 + df1 * x), df2 / 2.0, df1 / 2.0));
1371 }
1372
SCYTHE_DISTFUN_MATRIX(pf,double,SCYTHE_ARGSET (df1,df2),double df1,double df2)1373 SCYTHE_DISTFUN_MATRIX(pf, double, SCYTHE_ARGSET(df1, df2), double df1, double df2)
1374
1375 /* PDFs */
1376
1377 /*! \brief The F density function.
1378 *
1379 * Computes the value of the F probability density function with
1380 * \a df1 and \a df2 degrees of freedom, at the desired quantile \a
1381 * x.
1382 *
1383 * It is also possible to call this function with a Matrix of
1384 * doubles as its first argument. In this case the function will
1385 * return a Matrix of doubles of the same dimension as \a x,
1386 * containing the result of evaluating this function at each value
1387 * in \a x, given the remaining fixed parameters. By default, the
1388 * returned Matrix will be concrete and have the same matrix_order
1389 * as \a x, but you may invoke a generalized version of the function
1390 * with an explicit template call.
1391 *
1392 * \param x The desired quantile.
1393 * \param df1 The non-negative degrees of freedom for the
1394 * \f$\chi^2\f$ variate in the nominator of the F statistic.
1395 * \param df2 The non-negative degrees of freedom for the
1396 * \f$\chi^2\f$ variate in the denominator of the F statistic.
1397 *
1398 * \see df(double x, double df1, double df2)
1399 * \see rng::rf(double df1, double df2)
1400 *
1401 * \throw scythe_invalid_arg (Level 1)
1402 * \throw scythe_range_error (Level 1)
1403 * \throw scythe_precision_error (Level 1)
1404 * \throw scythe_convergence_error (Level 1)
1405 */
1406 inline double
1407 df(double x, double df1, double df2)
1408 {
1409 double dens;
1410
1411 SCYTHE_CHECK_10(df1 <= 0 || df2 <= 0, scythe_invalid_arg,
1412 "df1 or df2 <= 0");
1413
1414 if (x <= 0)
1415 return 0;
1416
1417 double f = 1 / (df2 + x * df1);
1418 double q = df2 * f;
1419 double p = x * df1 * f;
1420
1421 if (df1 >= 2) {
1422 f = df1 * q / 2;
1423 dens = dbinom_raw((df1 - 2) / 2,(df1 + df2 - 2) / 2, p, q);
1424 } else {
1425 f = (df1 * df1 * q) /(2 * p * (df1 + df2));
1426 dens = dbinom_raw(df1 / 2,(df1 + df2)/ 2, p, q);
1427 }
1428
1429 return f*dens;
1430 }
1431
SCYTHE_DISTFUN_MATRIX(df,double,SCYTHE_ARGSET (df1,df2),double df1,double df2)1432 SCYTHE_DISTFUN_MATRIX(df, double, SCYTHE_ARGSET(df1, df2), double df1, double df2)
1433
1434 /**** The Gamma Distribution ****/
1435
1436 /* CDFs */
1437 /*! \brief The gamma distribution function.
1438 *
1439 * Computes the value of the gamma cumulative distribution
1440 * function with given \a shape and \a scale, at the desired quantile
1441 * \a x.
1442 *
1443 * It is also possible to call this function with a Matrix of
1444 * doubles as its first argument. In this case the function will
1445 * return a Matrix of doubles of the same dimension as \a x,
1446 * containing the result of evaluating this function at each value
1447 * in \a x, given the remaining fixed parameters. By default, the
1448 * returned Matrix will be concrete and have the same matrix_order
1449 * as \a x, but you may invoke a generalized version of the function
1450 * with an explicit template call.
1451 *
1452 * \param x The desired quantile.
1453 * \param shape The non-negative shape of the distribution.
1454 * \param scale The non-negative scale of the distribution.
1455 *
1456 * \see dgamma(double x, double shape, double scale)
1457 * \see rng::rgamma(double shape, double scale)
1458 * \see gammafn(double x)
1459 * \see lngammafn(double x)
1460 *
1461 * \throw scythe_invalid_arg (Level 1)
1462 * \throw scythe_range_error (Level 1)
1463 * \throw scythe_precision_error (Level 1)
1464 * \throw scythe_convergence_error (Level 1)
1465 */
1466 inline double
1467 pgamma (double x, double shape, double scale)
1468 {
1469 const double xbig = 1.0e+8, xlarge = 1.0e+37,
1470 alphlimit = 1000.;/* normal approx. for shape > alphlimit */
1471
1472 int lower_tail = 1;
1473
1474 double pn1, pn2, pn3, pn4, pn5, pn6, arg, a, b, c, an, osum, sum;
1475 long n;
1476 int pearson;
1477
1478 /* check that we have valid values for x and shape */
1479
1480 SCYTHE_CHECK_10(shape <= 0. || scale <= 0., scythe_invalid_arg,
1481 "shape or scale <= 0");
1482
1483 x /= scale;
1484
1485 if (x <= 0.)
1486 return 0.0;
1487
1488 /* use a normal approximation if shape > alphlimit */
1489
1490 if (shape > alphlimit) {
1491 pn1 = std::sqrt(shape) * 3. * (std::pow(x/shape, 1./3.) + 1.
1492 / (9. * shape) - 1.);
1493 return pnorm(pn1, 0., 1.);
1494 }
1495
1496 /* if x is extremely large __compared to shape__ then return 1 */
1497
1498 if (x > xbig * shape)
1499 return 1.0;
1500
1501 if (x <= 1. || x < shape) {
1502 pearson = 1;/* use pearson's series expansion. */
1503 arg = shape * std::log(x) - x - lngammafn(shape + 1.);
1504 c = 1.;
1505 sum = 1.;
1506 a = shape;
1507 do {
1508 a += 1.;
1509 c *= x / a;
1510 sum += c;
1511 } while (c > DBL_EPSILON);
1512 arg += std::log(sum);
1513 }
1514 else { /* x >= max( 1, shape) */
1515 pearson = 0;/* use a continued fraction expansion */
1516 arg = shape * std::log(x) - x - lngammafn(shape);
1517 a = 1. - shape;
1518 b = a + x + 1.;
1519 pn1 = 1.;
1520 pn2 = x;
1521 pn3 = x + 1.;
1522 pn4 = x * b;
1523 sum = pn3 / pn4;
1524 for (n = 1; ; n++) {
1525 a += 1.;/* = n+1 -shape */
1526 b += 2.;/* = 2(n+1)-shape+x */
1527 an = a * n;
1528 pn5 = b * pn3 - an * pn1;
1529 pn6 = b * pn4 - an * pn2;
1530 if (std::fabs(pn6) > 0.) {
1531 osum = sum;
1532 sum = pn5 / pn6;
1533 if (std::fabs(osum - sum) <= DBL_EPSILON * std::min(1., sum))
1534 break;
1535 }
1536 pn1 = pn3;
1537 pn2 = pn4;
1538 pn3 = pn5;
1539 pn4 = pn6;
1540 if (std::fabs(pn5) >= xlarge) {
1541 /* re-scale terms in continued fraction if they are large */
1542 pn1 /= xlarge;
1543 pn2 /= xlarge;
1544 pn3 /= xlarge;
1545 pn4 /= xlarge;
1546 }
1547 }
1548 arg += std::log(sum);
1549 }
1550
1551 lower_tail = (lower_tail == pearson);
1552
1553 sum = std::exp(arg);
1554
1555 return (lower_tail) ? sum : 1 - sum;
1556 }
1557
SCYTHE_DISTFUN_MATRIX(pgamma,double,SCYTHE_ARGSET (shape,scale),double shape,double scale)1558 SCYTHE_DISTFUN_MATRIX(pgamma, double, SCYTHE_ARGSET(shape, scale), double shape, double scale)
1559
1560 /* PDFs */
1561 /*! \brief The gamma density function.
1562 *
1563 * Computes the value of the gamma probability density
1564 * function with given \a shape and \a scale, at the desired quantile
1565 * \a x.
1566 *
1567 * It is also possible to call this function with a Matrix of
1568 * doubles as its first argument. In this case the function will
1569 * return a Matrix of doubles of the same dimension as \a x,
1570 * containing the result of evaluating this function at each value
1571 * in \a x, given the remaining fixed parameters. By default, the
1572 * returned Matrix will be concrete and have the same matrix_order
1573 * as \a x, but you may invoke a generalized version of the function
1574 * with an explicit template call.
1575 *
1576 * \param x The desired quantile.
1577 * \param shape The non-negative shape of the distribution.
1578 * \param scale The non-negative scale of the distribution.
1579 *
1580 * \see pgamma(double x, double shape, double scale)
1581 * \see rng::rgamma(double shape, double scale)
1582 * \see gammafn(double x)
1583 * \see lngammafn(double x)
1584 *
1585 * \throw scythe_invalid_arg (Level 1)
1586 * \throw scythe_range_error (Level 1)
1587 * \throw scythe_precision_error (Level 1)
1588 * \throw scythe_convergence_error (Level 1)
1589 */
1590 inline double
1591 dgamma(double x, double shape, double scale)
1592 {
1593 SCYTHE_CHECK_10(shape <= 0 || scale <= 0,scythe_invalid_arg,
1594 "shape or scale <= 0");
1595
1596 if (x < 0)
1597 return 0.0;
1598
1599 if (x == 0) {
1600 SCYTHE_CHECK_10(shape < 1,scythe_invalid_arg,
1601 "x == 0 and shape < 1");
1602
1603 if (shape > 1)
1604 return 0.0;
1605
1606 return 1 / scale;
1607 }
1608
1609 if (shape < 1) {
1610 double pr = dpois_raw(shape, x/scale);
1611 return pr * shape / x;
1612 }
1613
1614 /* else shape >= 1 */
1615 double pr = dpois_raw(shape - 1, x / scale);
1616 return pr / scale;
1617 }
1618
SCYTHE_DISTFUN_MATRIX(dgamma,double,SCYTHE_ARGSET (shape,scale),double shape,double scale)1619 SCYTHE_DISTFUN_MATRIX(dgamma, double, SCYTHE_ARGSET(shape, scale), double shape, double scale)
1620
1621 /**** The Logistic Distribution ****/
1622
1623 /* CDFs */
1624 /*! \brief The logistic distribution function.
1625 *
1626 * Computes the value of the logistic cumulative distribution
1627 * function with given \a location and \a scale, at the desired
1628 * quantile \a x.
1629 *
1630 * It is also possible to call this function with a Matrix of
1631 * doubles as its first argument. In this case the function will
1632 * return a Matrix of doubles of the same dimension as \a x,
1633 * containing the result of evaluating this function at each value
1634 * in \a x, given the remaining fixed parameters. By default, the
1635 * returned Matrix will be concrete and have the same matrix_order
1636 * as \a x, but you may invoke a generalized version of the function
1637 * with an explicit template call.
1638 *
1639 * \param x The desired quantile.
1640 * \param location The location of the distribution.
1641 * \param scale The positive scale of the distribution.
1642 *
1643 * \see dlogis(double x, double location, double scale)
1644 * \see rng::rlogis(double location, double scale)
1645 *
1646 * \throw scythe_invalid_arg (Level 1)
1647 */
1648 inline double
1649 plogis (double x, double location, double scale)
1650 {
1651 SCYTHE_CHECK_10(scale <= 0.0, scythe_invalid_arg, "scale <= 0");
1652
1653 double X = (x-location) / scale;
1654
1655 X = std::exp(-X);
1656
1657 return 1 / (1+X);
1658 }
1659
SCYTHE_DISTFUN_MATRIX(plogis,double,SCYTHE_ARGSET (location,scale),double location,double scale)1660 SCYTHE_DISTFUN_MATRIX(plogis, double, SCYTHE_ARGSET(location, scale), double location, double scale)
1661
1662 /* PDFs */
1663 /*! \brief The logistic density function.
1664 *
1665 * Computes the value of the logistic probability density
1666 * function with given \a location and \a scale, at the desired
1667 * quantile \a x.
1668 *
1669 * It is also possible to call this function with a Matrix of
1670 * doubles as its first argument. In this case the function will
1671 * return a Matrix of doubles of the same dimension as \a x,
1672 * containing the result of evaluating this function at each value
1673 * in \a x, given the remaining fixed parameters. By default, the
1674 * returned Matrix will be concrete and have the same matrix_order
1675 * as \a x, but you may invoke a generalized version of the function
1676 * with an explicit template call.
1677 *
1678 * \param x The desired quantile.
1679 * \param location The location of the distribution.
1680 * \param scale The positive scale of the distribution.
1681 *
1682 * \see plogis(double x, double location, double scale)
1683 * \see rng::rlogis(double location, double scale)
1684 *
1685 * \throw scythe_invalid_arg (Level 1)
1686 */
1687 inline double
1688 dlogis(double x, double location, double scale)
1689 {
1690 SCYTHE_CHECK_10(scale <= 0.0, scythe_invalid_arg, "scale <= 0");
1691
1692 double X = (x - location) / scale;
1693 double e = std::exp(-X);
1694 double f = 1.0 + e;
1695
1696 return e / (scale * f * f);
1697 }
1698
SCYTHE_DISTFUN_MATRIX(dlogis,double,SCYTHE_ARGSET (location,scale),double location,double scale)1699 SCYTHE_DISTFUN_MATRIX(dlogis, double, SCYTHE_ARGSET(location, scale), double location, double scale)
1700
1701 /**** The Log Normal Distribution ****/
1702
1703 /* CDFs */
1704 /*! \brief The log-normal distribution function.
1705 *
1706 * Computes the value of the log-normal cumulative distribution
1707 * function with mean \a logmean and standard
1708 * deviation \a logsd, at the desired quantile \a x.
1709 *
1710 * It is also possible to call this function with a Matrix of
1711 * doubles as its first argument. In this case the function will
1712 * return a Matrix of doubles of the same dimension as \a x,
1713 * containing the result of evaluating this function at each value
1714 * in \a x, given the remaining fixed parameters. By default, the
1715 * returned Matrix will be concrete and have the same matrix_order
1716 * as \a x, but you may invoke a generalized version of the function
1717 * with an explicit template call.
1718 *
1719 * \param x The desired quantile.
1720 * \param logmean The mean of the distribution.
1721 * \param logsd The positive standard deviation of the distribution.
1722 *
1723 * \see dlnorm(double x, double logmean, double logsd)
1724 * \see rng::rlnorm(double logmean, double logsd)
1725 * \see pnorm(double x, double logmean, double logsd)
1726 *
1727 * \throw scythe_invalid_arg (Level 1)
1728 * \throw scythe_convergence_error (Level 1)
1729 */
1730 inline double
1731 plnorm (double x, double logmean, double logsd)
1732 {
1733 SCYTHE_CHECK_10(logsd <= 0, scythe_invalid_arg, "logsd <= 0");
1734
1735 if (x > 0)
1736 return pnorm(std::log(x), logmean, logsd);
1737
1738 return 0;
1739 }
1740
SCYTHE_DISTFUN_MATRIX(plnorm,double,SCYTHE_ARGSET (logmean,logsd),double logmean,double logsd)1741 SCYTHE_DISTFUN_MATRIX(plnorm, double, SCYTHE_ARGSET(logmean, logsd), double logmean, double logsd)
1742
1743 /* PDFs */
1744 /*! \brief The log-normal density function.
1745 *
1746 * Computes the value of the log-normal probability density
1747 * function with mean \a logmean and standard
1748 * deviation \a logsd, at the desired quantile \a x.
1749 *
1750 * It is also possible to call this function with a Matrix of
1751 * doubles as its first argument. In this case the function will
1752 * return a Matrix of doubles of the same dimension as \a x,
1753 * containing the result of evaluating this function at each value
1754 * in \a x, given the remaining fixed parameters. By default, the
1755 * returned Matrix will be concrete and have the same matrix_order
1756 * as \a x, but you may invoke a generalized version of the function
1757 * with an explicit template call.
1758 *
1759 * \param x The desired quantile.
1760 * \param logmean The mean of the distribution.
1761 * \param logsd The positive standard deviation of the distribution.
1762 *
1763 * \see plnorm(double x, double logmean, double logsd)
1764 * \see rng::rlnorm(double logmean, double logsd)
1765 * \see dnorm(double x, double logmean, double logsd)
1766 *
1767 * \throw scythe_invalid_arg (Level 1)
1768 */
1769 inline double
1770 dlnorm(double x, double logmean, double logsd)
1771 {
1772 SCYTHE_CHECK_10(logsd <= 0, scythe_invalid_arg, "logsd <= 0");
1773
1774 if (x == 0)
1775 return 0;
1776
1777 double y = (std::log(x) - logmean) / logsd;
1778
1779 return (1 / (std::sqrt(2 * M_PI))) * std::exp(-0.5 * y * y) / (x * logsd);
1780 }
1781
SCYTHE_DISTFUN_MATRIX(dlnorm,double,SCYTHE_ARGSET (logmean,logsd),double logmean,double logsd)1782 SCYTHE_DISTFUN_MATRIX(dlnorm, double, SCYTHE_ARGSET(logmean, logsd), double logmean, double logsd)
1783
1784 /**** The Negative Binomial Distribution ****/
1785
1786 /* CDFs */
1787 /*! \brief The negative binomial distribution function.
1788 *
1789 * Computes the value of the negative binomial cumulative distribution
1790 * function with \a n target number of successful trials and \a p
1791 * probability of success on each trial, at the desired quantile \a
1792 * x.
1793 *
1794 * It is also possible to call this function with a Matrix of
1795 * doubles as its first argument. In this case the function will
1796 * return a Matrix of doubles of the same dimension as \a x,
1797 * containing the result of evaluating this function at each value
1798 * in \a x, given the remaining fixed parameters. By default, the
1799 * returned Matrix will be concrete and have the same matrix_order
1800 * as \a x, but you may invoke a generalized version of the function
1801 * with an explicit template call.
1802 *
1803 * \param x The desired non-negative, integer, quantile.
1804 * \param n The positive target number of successful trials
1805 * (dispersion parameter).
1806 * \param p The probability of success on each trial.
1807 *
1808 * \see dnbinom(unsigned int x, double n, double p)
1809 * \see rng::rnbinom(double n, double p)
1810 *
1811 * \throw scythe_invalid_arg (Level 1)
1812 * \throw scythe_range_error (Level 1)
1813 * \throw scythe_precision_error (Level 1)
1814 */
1815 inline double
1816 pnbinom(unsigned int x, double n, double p)
1817 {
1818 SCYTHE_CHECK_10(n == 0 || p <= 0 || p >= 1, scythe_invalid_arg,
1819 "n == 0 or p not in (0,1)");
1820
1821 return pbeta(p, n, x + 1);
1822 }
1823
SCYTHE_DISTFUN_MATRIX(pnbinom,unsigned int,SCYTHE_ARGSET (n,p),double n,double p)1824 SCYTHE_DISTFUN_MATRIX(pnbinom, unsigned int, SCYTHE_ARGSET(n, p), double n, double p)
1825
1826 /* PDFs */
1827 /*! \brief The negative binomial density function.
1828 *
1829 * Computes the value of the negative binomial probability density
1830 * function with \a n target number of successful trials and \a p
1831 * probability of success on each trial, at the desired quantile \a
1832 * x.
1833 *
1834 * It is also possible to call this function with a Matrix of
1835 * doubles as its first argument. In this case the function will
1836 * return a Matrix of doubles of the same dimension as \a x,
1837 * containing the result of evaluating this function at each value
1838 * in \a x, given the remaining fixed parameters. By default, the
1839 * returned Matrix will be concrete and have the same matrix_order
1840 * as \a x, but you may invoke a generalized version of the function
1841 * with an explicit template call.
1842 *
1843 * \param x The desired non-negative, integer, quantile.
1844 * \param n The positive target number of successful trials
1845 * (dispersion parameter).
1846 * \param p The probability of success on each trial.
1847 *
1848 * \see dnbinom(unsigned int x, double n, double p)
1849 * \see rng::rnbinom(double n, double p)
1850 *
1851 * \throw scythe_invalid_arg (Level 1)
1852 * \throw scythe_range_error (Level 1)
1853 * \throw scythe_precision_error (Level 1)
1854 */
1855 inline double
1856 dnbinom(unsigned int x, double n, double p)
1857 {
1858 SCYTHE_CHECK_10(n == 0 || p <= 0 || p >= 1, scythe_invalid_arg,
1859 "n == 0 or p not in (0,1)");
1860
1861 double prob = dbinom_raw(n, x + n, p, 1 - p);
1862 double P = (double) n / (n + x);
1863
1864 return P * prob;
1865 }
1866
SCYTHE_DISTFUN_MATRIX(dnbinom,unsigned int,SCYTHE_ARGSET (n,p),double n,double p)1867 SCYTHE_DISTFUN_MATRIX(dnbinom, unsigned int, SCYTHE_ARGSET(n, p), double n, double p)
1868
1869 /**** The Normal Distribution ****/
1870
1871 /* CDFs */
1872 /*! \brief The normal distribution function.
1873 *
1874 * Computes the value of the normal cumulative distribution
1875 * function with given \a mean and standard deviation \a sd, at the
1876 * desired quantile \a x.
1877 *
1878 * It is also possible to call this function with a Matrix of
1879 * doubles as its first argument. In this case the function will
1880 * return a Matrix of doubles of the same dimension as \a x,
1881 * containing the result of evaluating this function at each value
1882 * in \a x, given the remaining fixed parameters. By default, the
1883 * returned Matrix will be concrete and have the same matrix_order
1884 * as \a x, but you may invoke a generalized version of the function
1885 * with an explicit template call.
1886 *
1887 * \param x The desired quantile.
1888 * \param mean The mean of the distribution.
1889 * \param sd The positive standard deviation of the distribution.
1890 *
1891 * \see dnorm(double x, double mean, double sd)
1892 * \see rng::rnorm(double mean, double sd)
1893 *
1894 * \throw scythe_invalid_arg (Level 1)
1895 * \throw scythe_convergence_error (Level 1)
1896 */
1897 inline double
1898 pnorm (double x, double mean, double sd)
1899
1900 {
1901 SCYTHE_CHECK_10(sd <= 0, scythe_invalid_arg,
1902 "negative standard deviation");
1903
1904 return pnorm1((x - mean) / sd, true, false);
1905 }
1906
SCYTHE_DISTFUN_MATRIX(pnorm,double,SCYTHE_ARGSET (mean,sd),double mean,double sd)1907 SCYTHE_DISTFUN_MATRIX(pnorm, double, SCYTHE_ARGSET(mean, sd), double mean, double sd)
1908
1909
1910 /* PDFs */
1911 /*! \brief The normal density function.
1912 *
1913 * Computes the value of the normal probability density
1914 * function with given \a mean and standard deviation \a sd, at the
1915 * desired quantile \a x.
1916 *
1917 * It is also possible to call this function with a Matrix of
1918 * doubles as its first argument. In this case the function will
1919 * return a Matrix of doubles of the same dimension as \a x,
1920 * containing the result of evaluating this function at each value
1921 * in \a x, given the remaining fixed parameters. By default, the
1922 * returned Matrix will be concrete and have the same matrix_order
1923 * as \a x, but you may invoke a generalized version of the function
1924 * with an explicit template call.
1925 *
1926 * \param x The desired quantile.
1927 * \param mean The mean of the distribution.
1928 * \param sd The positive standard deviation of the distribution.
1929 *
1930 * \see pnorm(double x, double mean, double sd)
1931 * \see rng::rnorm(double mean, double sd)
1932 *
1933 * \throw scythe_invalid_arg (Level 1)
1934 */
1935 inline double
1936 dnorm(double x, double mean, double sd)
1937 {
1938 SCYTHE_CHECK_10(sd <= 0, scythe_invalid_arg,
1939 "negative standard deviation");
1940
1941 double X = (x - mean) / sd;
1942
1943 return (M_1_SQRT_2PI * std::exp(-0.5 * X * X) / sd);
1944 }
1945
SCYTHE_DISTFUN_MATRIX(dnorm,double,SCYTHE_ARGSET (mean,sd),double mean,double sd)1946 SCYTHE_DISTFUN_MATRIX(dnorm, double, SCYTHE_ARGSET(mean, sd), double mean, double sd)
1947
1948
1949 /* Return the natural log of the normal PDF */
1950 /*! \brief The natural log of normal density function.
1951 *
1952 * Computes the value of the natural log of the normal probability
1953 * density function with given \a mean and standard deviation \a sd,
1954 * at the desired quantile \a x.
1955 *
1956 * It is also possible to call this function with a Matrix of
1957 * doubles as its first argument. In this case the function will
1958 * return a Matrix of doubles of the same dimension as \a x,
1959 * containing the result of evaluating this function at each value
1960 * in \a x, given the remaining fixed parameters. By default, the
1961 * returned Matrix will be concrete and have the same matrix_order
1962 * as \a x, but you may invoke a generalized version of the function
1963 * with an explicit template call.
1964 *
1965 * \param x The desired quantile.
1966 * \param mean The mean of the distribution.
1967 * \param sd The positive standard deviation of the distribution.
1968 *
1969 * \see dnorm(double x, double mean, double sd)
1970 * \see pnorm(double x, double mean, double sd)
1971 * \see rng::rnorm(double mean, double sd)
1972 *
1973 * \throw scythe_invalid_arg (Level 1)
1974 */
1975 inline double
1976 lndnorm (double x, double mean, double sd)
1977 {
1978 SCYTHE_CHECK_10(sd <= 0, scythe_invalid_arg,
1979 "negative standard deviation");
1980
1981 double X = (x - mean) / sd;
1982
1983 return -(M_LN_SQRT_2PI + 0.5 * X * X + std::log(sd));
1984 }
1985
SCYTHE_DISTFUN_MATRIX(lndnorm,double,SCYTHE_ARGSET (mean,sd),double mean,double sd)1986 SCYTHE_DISTFUN_MATRIX(lndnorm, double, SCYTHE_ARGSET(mean, sd), double mean, double sd)
1987
1988 /* Quantile functions */
1989 /*! \brief The standard normal quantile function.
1990 *
1991 * Computes the value of the standard normal quantile function
1992 * at the desired probability \a in_p.
1993 *
1994 * It is also possible to call this function with a Matrix of
1995 * doubles as its first argument. In this case the function will
1996 * return a Matrix of doubles of the same dimension as \a x,
1997 * containing the result of evaluating this function at each value
1998 * in \a x, given the remaining fixed parameters. By default, the
1999 * returned Matrix will be concrete and have the same matrix_order
2000 * as \a x, but you may invoke a generalized version of the function
2001 * with an explicit template call.
2002 *
2003 * \param in_p The desired probability.
2004 *
2005 * \see pnorm(double x, double mean, double sd)
2006 * \see dnorm(double x, double mean, double sd)
2007 * \see rng::rnorm(double mean, double sd)
2008 *
2009 * \throw scythe_invalid_arg (Level 1)
2010 */
2011 inline double
2012 qnorm1 (double in_p)
2013 {
2014 double p0 = -0.322232431088;
2015 double q0 = 0.0993484626060;
2016 double p1 = -1.0;
2017 double q1 = 0.588581570495;
2018 double p2 = -0.342242088547;
2019 double q2 = 0.531103462366;
2020 double p3 = -0.0204231210245;
2021 double q3 = 0.103537752850;
2022 double p4 = -0.453642210148e-4;
2023 double q4 = 0.38560700634e-2;
2024 double xp = 0.0;
2025 double p = in_p;
2026
2027 if (p > 0.5)
2028 p = 1 - p;
2029
2030 SCYTHE_CHECK_10(p < 10e-20, scythe_range_error,
2031 "p outside accuracy limit");
2032
2033 if (p == 0.5)
2034 return xp;
2035
2036 double y = std::sqrt (std::log (1.0 / std::pow (p, 2)));
2037 xp = y + ((((y * p4 + p3) * y + p2) * y + p1) * y + p0) /
2038 ((((y * q4 + q3) * y + q2) * y + q1) * y + q0);
2039
2040 if (in_p < 0.5)
2041 xp = -1 * xp;
2042
2043 return xp;
2044 }
2045
SCYTHE_DISTFUN_MATRIX(qnorm1,double,in_p,double in_p)2046 SCYTHE_DISTFUN_MATRIX(qnorm1, double, in_p, double in_p)
2047
2048 /**** The Poisson Distribution ****/
2049
2050 /* CDFs */
2051 /*! \brief The Poisson distribution function.
2052 *
2053 * Computes the value of the Poisson cumulative distribution
2054 * function with expected number of occurrences \a lambda, at the
2055 * desired quantile \a x.
2056 *
2057 * It is also possible to call this function with a Matrix of
2058 * doubles as its first argument. In this case the function will
2059 * return a Matrix of doubles of the same dimension as \a x,
2060 * containing the result of evaluating this function at each value
2061 * in \a x, given the remaining fixed parameters. By default, the
2062 * returned Matrix will be concrete and have the same matrix_order
2063 * as \a x, but you may invoke a generalized version of the function
2064 * with an explicit template call.
2065 *
2066 * \param x The desired integer quantile.
2067 * \param lambda The expected positive number of occurrences.
2068 *
2069 * \see dpois(unsigned int x, double lambda)
2070 * \see rng::rpois(double lambda)
2071 *
2072 * \throws scythe_invalid_arg (Level 1)
2073 * \throw scythe_range_error (Level 1)
2074 * \throw scythe_precision_error (Level 1)
2075 * \throw scythe_convergence_error (Level 1)
2076 */
2077 inline double
2078 ppois(unsigned int x, double lambda)
2079 {
2080 SCYTHE_CHECK_10(lambda<=0.0, scythe_invalid_arg, "lambda <= 0");
2081
2082 if (lambda == 1)
2083 return 1;
2084
2085 return 1 - pgamma(lambda, x + 1, 1.0);
2086 }
2087
SCYTHE_DISTFUN_MATRIX(ppois,unsigned int,lambda,double lambda)2088 SCYTHE_DISTFUN_MATRIX(ppois, unsigned int, lambda, double lambda)
2089
2090 /* PDFs */
2091 /*! \brief The Poisson density function.
2092 *
2093 * Computes the value of the Poisson probability density
2094 * function with expected number of occurrences \a lambda, at the
2095 * desired quantile \a x.
2096 *
2097 * It is also possible to call this function with a Matrix of
2098 * doubles as its first argument. In this case the function will
2099 * return a Matrix of doubles of the same dimension as \a x,
2100 * containing the result of evaluating this function at each value
2101 * in \a x, given the remaining fixed parameters. By default, the
2102 * returned Matrix will be concrete and have the same matrix_order
2103 * as \a x, but you may invoke a generalized version of the function
2104 * with an explicit template call.
2105 *
2106 * \param x The desired integer quantile.
2107 * \param lambda The expected positive number of occurrences.
2108 *
2109 * \see ppois(unsigned int x, double lambda)
2110 * \see rng::rpois(double lambda)
2111 *
2112 * \throws scythe_invalid_arg (Level 1)
2113 */
2114 inline double
2115 dpois(unsigned int x, double lambda)
2116 {
2117 SCYTHE_CHECK_10(lambda<=0.0, scythe_invalid_arg, "lambda <= 0");
2118
2119 // compute log(x!)
2120 double xx = x+1;
2121 double cof[6] = {
2122 76.18009172947146, -86.50532032941677,
2123 24.01409824083091, -1.231739572450155,
2124 0.1208650973866179e-2, -0.5395239384953e-5
2125 };
2126 double y = xx;
2127 double tmp = xx + 5.5 - (xx + 0.5) * std::log(xx + 5.5);
2128 double ser = 1.000000000190015;
2129 for (int j = 0; j <= 5; j++) {
2130 ser += (cof[j] / ++y);
2131 }
2132 double lnfactx = std::log(2.5066282746310005 * ser / xx) - tmp;
2133
2134 return (std::exp( -1*lnfactx + x * std::log(lambda) - lambda));
2135 }
2136
SCYTHE_DISTFUN_MATRIX(dpois,unsigned int,lambda,double lambda)2137 SCYTHE_DISTFUN_MATRIX(dpois, unsigned int, lambda, double lambda)
2138
2139 /**** The t Distribution ****/
2140
2141 /* CDFs */
2142 /*! \brief The Student t distribution function.
2143 *
2144 * Computes the value of the Student t cumulative distribution
2145 * function with \a n degrees of freedom, at the desired quantile
2146 * \a x.
2147 *
2148 * It is also possible to call this function with a Matrix of
2149 * doubles as its first argument. In this case the function will
2150 * return a Matrix of doubles of the same dimension as \a x,
2151 * containing the result of evaluating this function at each value
2152 * in \a x, given the remaining fixed parameters. By default, the
2153 * returned Matrix will be concrete and have the same matrix_order
2154 * as \a x, but you may invoke a generalized version of the function
2155 * with an explicit template call.
2156 *
2157 * \param x The desired quantile.
2158 * \param n The positive degrees of freedom of the distribution.
2159 *
2160 * \see dt(double x, bool b1, bool b2)
2161 * \see rng::rt1(double mu, double sigma2, double nu)
2162 *
2163 * \throw scythe_invalid_arg (Level 1)
2164 * \throw scythe_convergence_error (Level 1)
2165 * \throw scythe_range_error (Level 1)
2166 * \throw scythe_precision_error (Level 1)
2167 */
2168 inline double
2169 pt(double x, double n)
2170 {
2171 double val;
2172
2173 SCYTHE_CHECK_10(n <= 0, scythe_invalid_arg, "n <= 0");
2174
2175 if (n > 4e5) {
2176 val = 1/(4*n);
2177 return pnorm1(x * (1 - val) / ::sqrt(1 + x * x * 2. * val),
2178 true, false);
2179 }
2180
2181 val = pbeta(n / (n + x * x), n / 2.0, 0.5);
2182
2183 val /= 2;
2184
2185 if (x <= 0)
2186 return val;
2187 else
2188 return 1 - val;
2189 }
2190
SCYTHE_DISTFUN_MATRIX(pt,double,n,double n)2191 SCYTHE_DISTFUN_MATRIX(pt, double, n, double n)
2192
2193 /* PDFs */
2194 /*! \brief The Student t distribution function.
2195 *
2196 * Computes the value of the Student t cumulative distribution
2197 * function with \a n degrees of freedom, at the desired quantile
2198 * \a x.
2199 *
2200 * It is also possible to call this function with a Matrix of
2201 * doubles as its first argument. In this case the function will
2202 * return a Matrix of doubles of the same dimension as \a x,
2203 * containing the result of evaluating this function at each value
2204 * in \a x, given the remaining fixed parameters. By default, the
2205 * returned Matrix will be concrete and have the same matrix_order
2206 * as \a x, but you may invoke a generalized version of the function
2207 * with an explicit template call.
2208 *
2209 * \param x The desired quantile.
2210 * \param n The positive degrees of freedom of the distribution.
2211 *
2212 * \see pt(double x, bool b1, bool b2)
2213 * \see rng::rt1(double mu, double sigma2, double nu)
2214 *
2215 * \throw scythe_invalid_arg (Level 1)
2216 * \throw scythe_range_error (Level 1)
2217 * \throw scythe_precision_error (Level 1)
2218 */
2219 inline double
2220 dt(double x, double n)
2221 {
2222 double u;
2223
2224 SCYTHE_CHECK_10(n <= 0, scythe_invalid_arg, "n <= 0");
2225
2226 double t = -bd0(n/2., (n + 1) / 2.)
2227 + stirlerr((n + 1) / 2.)
2228 - stirlerr(n / 2.);
2229 if(x*x > 0.2*n)
2230 u = std::log(1+x*x/n)*n/2;
2231 else
2232 u = -bd0(n/2., (n+x*x)/2.) + x*x/2;
2233
2234 return std::exp(t-u)/std::sqrt(2*M_PI*(1+x*x/n));
2235 }
2236
SCYTHE_DISTFUN_MATRIX(dt,double,n,double n)2237 SCYTHE_DISTFUN_MATRIX(dt, double, n, double n)
2238
2239 /* Returns the univariate Student-t density evaluated at x
2240 * with mean mu, scale sigma^2, and nu degrees of freedom.
2241 *
2242 * TODO: Do we want a pt for this distribution?
2243 */
2244
2245 /*! \brief The univariate Student t density function.
2246 *
2247 * Computes the value of the univariate Student t probability
2248 * density function with mean \a mu, variance \a sigma2,
2249 * and degrees of freedom \a nu, at the desired quantile \a x.
2250 *
2251 * It is also possible to call this function with a Matrix of
2252 * doubles as its first argument. In this case the function will
2253 * return a Matrix of doubles of the same dimension as \a x,
2254 * containing the result of evaluating this function at each value
2255 * in \a x, given the remaining fixed parameters. By default, the
2256 * returned Matrix will be concrete and have the same matrix_order
2257 * as \a x, but you may invoke a generalized version of the function
2258 * with an explicit template call.
2259 *
2260 * \param x The desired quantile.
2261 * \param mu The mean of the distribution.
2262 * \param sigma2 The variance of the distribution.
2263 * \param nu The degrees of freedom of the distribution.
2264 *
2265 * \see rng::rt1(double mu, double sigma2, double nu)
2266 * \see dt(double x, bool b1, bool b2)
2267 * \see pt(double x, bool b1, bool b2)
2268 *
2269 * \throw scythe_invalid_arg (Level 1)
2270 * \throw scythe_range_error (Level 1)
2271 * \throw scythe_precision_error (Level 1)
2272 */
2273 inline double
2274 dt1(double x, double mu, double sigma2, double nu)
2275 {
2276 double logdens = lngammafn((nu + 1.0) /2.0)
2277 - std::log(std::sqrt(nu * M_PI))
2278 - lngammafn(nu / 2.0) - std::log(std::sqrt(sigma2))
2279 - (nu + 1.0) / 2.0 * std::log(1.0 + (std::pow((x - mu), 2.0))
2280 / (nu * sigma2));
2281
2282 return(std::exp(logdens));
2283 }
2284
SCYTHE_DISTFUN_MATRIX(dt1,double,SCYTHE_ARGSET (mu,sigma2,nu),double mu,double sigma2,double nu)2285 SCYTHE_DISTFUN_MATRIX(dt1, double, SCYTHE_ARGSET(mu, sigma2, nu), double mu, double sigma2, double nu)
2286
2287 /* Returns the natural log of the univariate Student-t density
2288 * evaluated at x with mean mu, scale sigma^2, and nu
2289 * degrees of freedom
2290 */
2291
2292 /*! \brief The natural log of the univariate Student t density
2293 * function.
2294 *
2295 * Computes the value of the natural log of the univariate Student t
2296 * probability density function with mean \a mu, variance \a sigma2,
2297 * and degrees of freedom \a nu, at the desired quantile \a x.
2298 *
2299 * It is also possible to call this function with a Matrix of
2300 * doubles as its first argument. In this case the function will
2301 * return a Matrix of doubles of the same dimension as \a x,
2302 * containing the result of evaluating this function at each value
2303 * in \a x, given the remaining fixed parameters. By default, the
2304 * returned Matrix will be concrete and have the same matrix_order
2305 * as \a x, but you may invoke a generalized version of the function
2306 * with an explicit template call.
2307 *
2308 * \param x The desired quantile.
2309 * \param mu The mean of the distribution.
2310 * \param sigma2 The variance of the distribution.
2311 * \param nu The degrees of freedom of the distribution.
2312 *
2313 * \see rng::rt1(double mu, double sigma2, double nu)
2314 * \see dt(double x, bool b1, bool b2)
2315 * \see pt(double x, bool b1, bool b2)
2316 *
2317 * \throw scythe_invalid_arg (Level 1)
2318 * \throw scythe_range_error (Level 1)
2319 * \throw scythe_precision_error (Level 1)
2320 */
2321 inline double
2322 lndt1(double x, double mu, double sigma2, double nu)
2323 {
2324 double logdens = lngammafn((nu+1.0)/2.0)
2325 - std::log(std::sqrt(nu*M_PI))
2326 - lngammafn(nu/2.0) - std::log(std::sqrt(sigma2))
2327 - (nu+1.0)/2.0 * std::log(1.0 + (std::pow((x-mu),2.0))
2328 /(nu * sigma2));
2329
2330 return(logdens);
2331 }
2332
SCYTHE_DISTFUN_MATRIX(lndt1,double,SCYTHE_ARGSET (mu,sigma2,nu),double mu,double sigma2,double nu)2333 SCYTHE_DISTFUN_MATRIX(lndt1, double, SCYTHE_ARGSET(mu, sigma2, nu), double mu, double sigma2, double nu)
2334
2335 /**** The Uniform Distribution ****/
2336
2337 /* CDFs */
2338 /*! \brief The uniform distribution function.
2339 *
2340 * Computes the value of the uniform cumulative distribution
2341 * function evaluated on the interval [\a a, \a b], at the desired
2342 * quantile \a x.
2343 *
2344 * It is also possible to call this function with a Matrix of
2345 * doubles as its first argument. In this case the function will
2346 * return a Matrix of doubles of the same dimension as \a x,
2347 * containing the result of evaluating this function at each value
2348 * in \a x, given the remaining fixed parameters. By default, the
2349 * returned Matrix will be concrete and have the same matrix_order
2350 * as \a x, but you may invoke a generalized version of the function
2351 * with an explicit template call.
2352 *
2353 * \param x The desired quantile x.
2354 * \param a The lower end-point of the distribution.
2355 * \param b The upper end-point of the distribution.
2356 *
2357 * \see dunif(double x, double a, double b)
2358 * \see rng::runif()
2359 *
2360 * \throw scythe_invalid_arg (Level 1)
2361 */
2362 inline double
2363 punif(double x, double a, double b)
2364 {
2365 SCYTHE_CHECK_10(b <= a, scythe_invalid_arg, "b <= a");
2366
2367 if (x <= a)
2368 return 0.0;
2369
2370 if (x >= b)
2371 return 1.0;
2372
2373 return (x - a) / (b - a);
2374 }
2375
SCYTHE_DISTFUN_MATRIX(punif,double,SCYTHE_ARGSET (a,b),double a,double b)2376 SCYTHE_DISTFUN_MATRIX(punif, double, SCYTHE_ARGSET(a, b), double a, double b)
2377
2378 /* PDFs */
2379 /*! \brief The uniform density function.
2380 *
2381 * Computes the value of the uniform probability density
2382 * function evaluated on the interval [\a a, \a b], at the desired
2383 * quantile \a x.
2384 *
2385 * It is also possible to call this function with a Matrix of
2386 * doubles as its first argument. In this case the function will
2387 * return a Matrix of doubles of the same dimension as \a x,
2388 * containing the result of evaluating this function at each value
2389 * in \a x, given the remaining fixed parameters. By default, the
2390 * returned Matrix will be concrete and have the same matrix_order
2391 * as \a x, but you may invoke a generalized version of the function
2392 * with an explicit template call.
2393 *
2394 * \param x The desired quantile x.
2395 * \param a The lower end-point of the distribution.
2396 * \param b The upper end-point of the distribution.
2397 *
2398 * \see punif(double x, double a, double b)
2399 * \see rng::runif()
2400 *
2401 * \throw scythe_invalid_arg (Level 1)
2402 */
2403 inline double
2404 dunif(double x, double a, double b)
2405 {
2406 SCYTHE_CHECK_10(b <= a, scythe_invalid_arg, "b <= a");
2407
2408 if (a <= x && x <= b)
2409 return 1.0 / (b - a);
2410
2411 return 0.0;
2412 }
2413
SCYTHE_DISTFUN_MATRIX(dunif,double,SCYTHE_ARGSET (a,b),double a,double b)2414 SCYTHE_DISTFUN_MATRIX(dunif, double, SCYTHE_ARGSET(a, b), double a, double b)
2415
2416 /**** The Weibull Distribution ****/
2417
2418 /* CDFs */
2419 /*! \brief The Weibull distribution function.
2420 *
2421 * Computes the value of the Weibull cumulative distribution
2422 * function with given \a shape and \a scale, at the desired
2423 * quantile \a x.
2424 *
2425 * It is also possible to call this function with a Matrix of
2426 * doubles as its first argument. In this case the function will
2427 * return a Matrix of doubles of the same dimension as \a x,
2428 * containing the result of evaluating this function at each value
2429 * in \a x, given the remaining fixed parameters. By default, the
2430 * returned Matrix will be concrete and have the same matrix_order
2431 * as \a x, but you may invoke a generalized version of the function
2432 * with an explicit template call.
2433 *
2434 * \param x The desired quantile.
2435 * \param shape The positive shape of the distribution.
2436 * \param scale The positive scale of the distribution.
2437 *
2438 * \see dweibull(double x, double shape, double scale)
2439 * \see rng::rweibull(double shape, double scale)
2440 *
2441 * \throw scythe_invalid_arg (Level 1)
2442 */
2443 inline double
2444 pweibull(double x, double shape, double scale)
2445 {
2446 SCYTHE_CHECK_10(shape <= 0 || scale <= 0, scythe_invalid_arg,
2447 "shape or scale <= 0");
2448
2449 if (x <= 0)
2450 return 0.0;
2451
2452 return 1 - std::exp(-std::pow(x / scale, shape));
2453 }
2454
SCYTHE_DISTFUN_MATRIX(pweibull,double,SCYTHE_ARGSET (shape,scale),double shape,double scale)2455 SCYTHE_DISTFUN_MATRIX(pweibull, double, SCYTHE_ARGSET(shape, scale), double shape, double scale)
2456
2457 /* PDFs */
2458 /*! \brief The Weibull density function.
2459 *
2460 * Computes the value of the Weibull probability density
2461 * function with given \a shape and \a scale, at the desired
2462 * quantile \a x.
2463 *
2464 * It is also possible to call this function with a Matrix of
2465 * doubles as its first argument. In this case the function will
2466 * return a Matrix of doubles of the same dimension as \a x,
2467 * containing the result of evaluating this function at each value
2468 * in \a x, given the remaining fixed parameters. By default, the
2469 * returned Matrix will be concrete and have the same matrix_order
2470 * as \a x, but you may invoke a generalized version of the function
2471 * with an explicit template call.
2472 *
2473 * \param x The desired quantile.
2474 * \param shape The positive shape of the distribution.
2475 * \param scale The positive scale of the distribution.
2476 *
2477 * \see pweibull(double x, double shape, double scale)
2478 * \see rng::rweibull(double shape, double scale)
2479 *
2480 * \throw scythe_invalid_arg (Level 1)
2481 */
2482 inline double
2483 dweibull(double x, double shape, double scale)
2484 {
2485 SCYTHE_CHECK_10(shape <= 0 || scale <= 0, scythe_invalid_arg,
2486 "shape or scale <= 0");
2487
2488 if (x < 0)
2489 return 0.;
2490
2491 double tmp1 = std::pow(x / scale, shape - 1);
2492 double tmp2 = tmp1*(x / scale);
2493
2494 return shape * tmp1 * std::exp(-tmp2) / scale;
2495 }
2496
SCYTHE_DISTFUN_MATRIX(dweibull,double,SCYTHE_ARGSET (shape,scale),double shape,double scale)2497 SCYTHE_DISTFUN_MATRIX(dweibull, double, SCYTHE_ARGSET(shape, scale), double shape, double scale)
2498
2499 /* Multivariate Normal */
2500
2501 // TODO: distribution function. Plain old (non-logged) dmvnorm.
2502
2503 /*! \brief The natural log of the multivariate normal density
2504 * function.
2505 *
2506 * Computes the value of the natural log of the multivariate normal
2507 * probability density function with vector of mean \a mu and
2508 * variance-covariance matrix \a Sigma, at the vector of desired
2509 * quantiles \a x.
2510 *
2511 * \param x The vector of desired quantiles.
2512 * \param mu The vector of means.
2513 * \param Sigma The variance-covariance matrix.
2514 *
2515 * \see rng:rmvnorm(const Matrix<double, PO1, PS1>& mu, const Matrix<double, PO2, PS2>& sigma)
2516 *
2517 * \throw scythe_dimension_error (Level 1)
2518 * \throw scythe_conformation_error (Level 1)
2519 * \throw scythe_null_error (Level 1)
2520 */
2521 template <matrix_order O1, matrix_style S1,
2522 matrix_order O2, matrix_style S2,
2523 matrix_order O3, matrix_style S3>
2524 double lndmvn (const Matrix<double, O1, S1>& x,
2525 const Matrix<double, O2, S2>& mu,
2526 const Matrix<double, O3, S3>& Sigma)
2527 {
2528 SCYTHE_CHECK_10(! x.isColVector(), scythe_dimension_error,
2529 "x is not a column vector");
2530 SCYTHE_CHECK_10(! mu.isColVector(), scythe_dimension_error,
2531 "mu is not a column vector");
2532 SCYTHE_CHECK_10(! Sigma.isSquare(), scythe_dimension_error,
2533 "Sigma is not square");
2534 SCYTHE_CHECK_10(mu.rows()!=Sigma.rows() || x.rows()!=Sigma.rows(),
2535 scythe_conformation_error,
2536 "mu, x and Sigma have mismatched row lengths")
2537 int k = (int) mu.rows();
2538 return ( (-k/2.0)*std::log(2*M_PI) -0.5 * std::log(det(Sigma))
2539 -0.5 * (t(x - mu)) * invpd(Sigma) * (x-mu) )[0];
2540 }
2541
2542 } // end namespace scythe
2543
2544
2545 #endif /* SCYTHE_DISTRIBUTIONS_H */
2546