1 #include <config.h>
2 #include "DNT.h"
3
4 #include <cmath>
5
6 #include <JRmath.h>
7
8 using std::vector;
9
MU(vector<double const * > const & par)10 static inline double MU(vector<double const*> const &par) { return *par[0]; }
TAU(vector<double const * > const & par)11 static inline double TAU(vector<double const*> const &par) { return *par[1]; }
DF(vector<double const * > const & par)12 static inline double DF(vector<double const*> const &par) { return *par[2]; }
13
SIGMA(vector<double const * > const & par)14 static inline double SIGMA(vector<double const*> const &par) {
15 //Standard deviation of numerator
16 return 1/sqrt(TAU(par));
17 }
DELTA(vector<double const * > const & par)18 static inline double DELTA(vector<double const*> const &par) {
19 //Non-centrality parameter of scaled x
20 return MU(par)/SIGMA(par);
21 }
22
23 namespace jags {
24 namespace bugs {
25
DNT()26 DNT::DNT()
27 : RScalarDist("dnt", 3, DIST_UNBOUNDED)
28 {}
29
checkParameterValue(vector<double const * > const & par) const30 bool DNT::checkParameterValue (vector<double const *> const &par) const
31 {
32 return (TAU(par) > 0 && DF(par) > 0 && fabs(DELTA(par)) <= 37.62);
33 }
34
d(double x,PDFType type,vector<double const * > const & par,bool give_log) const35 double DNT::d(double x, PDFType type,
36 vector<double const *> const &par, bool give_log) const
37 {
38 x /= SIGMA(par);
39 if (give_log) {
40 return dnt(x, DF(par), DELTA(par), 1) - log(SIGMA(par));
41 }
42 else {
43 return dnt(x, DF(par), DELTA(par), 0) / SIGMA(par);
44 }
45 }
46
p(double x,vector<double const * > const & par,bool lower,bool use_log) const47 double DNT::p(double x, vector<double const *> const &par, bool lower,
48 bool use_log) const
49 {
50 return pnt(x/SIGMA(par), DF(par), DELTA(par), lower, use_log);
51 }
52
q(double p,vector<double const * > const & par,bool lower,bool log_p) const53 double DNT::q(double p, vector<double const *> const &par, bool lower,
54 bool log_p) const
55 {
56 return qnt(p, DF(par), DELTA(par), lower, log_p) * SIGMA(par);
57 }
58
r(vector<double const * > const & par,RNG * rng) const59 double DNT::r(vector<double const *> const &par, RNG *rng) const
60 {
61 double k = DF(par);
62 return rnorm(MU(par), SIGMA(par), rng)/sqrt(rchisq(k, rng)/k);
63 }
64 }
65 }
66