7 #include "boost/math/special_functions.hpp" 18 constexpr T inv_sqrt_2pi = T(0.3989422804014327);
22 constexpr T sqrt_two_over_pi(0.797884560802865);
26 constexpr T log_two_pi (1.83787706640935);
30 constexpr T log_two_over_pi (-0.451582705289455);
34 constexpr T log_pi (1.1447298858494);
46 template<
typename float_t>
47 float_t twiceFisher(float_t phi)
49 if ( (phi <= -1.0) || (phi >= 1.0) )
50 throw std::invalid_argument(
"error: phi was not between -1 and 1" );
52 return std::log(1.0 + phi) - std::log(1.0 - phi);
62 template<
typename float_t>
63 float_t invTwiceFisher(float_t psi)
65 float_t ans = (1.0 - std::exp(psi)) / ( -1.0 - std::exp(psi) );
67 if ( (ans <= -1.0) || (ans >= 1.0) )
68 std::cerr <<
"error: there was probably overflow for exp(psi) \n";
79 template<
typename float_t>
80 float_t logit(float_t p)
82 if ( (p <= 0.0) || (p >= 1.0))
83 std::cerr <<
"error: p was not between 0 and 1 \n";
85 return std::log(p) - std::log(1.0 - p);
94 template<
typename float_t>
95 float_t inv_logit(float_t r)
97 float_t ans = 1.0/( 1.0 + std::exp(-r) );
99 if ( (ans <= 0.0) || (ans >= 1.0))
100 std::cerr <<
"error: there was probably underflow for exp(-r) \n";
111 template<
typename float_t>
112 float_t log_inv_logit(float_t r)
114 if(r < -750.00 || r > 750.00) std::cerr <<
"warning: log_inv_logit might be under/over-flowing\n";
115 return -std::log(1.0 + std::exp(-r));
125 template<
typename float_t>
126 float_t log_sum_exp(float_t a, float_t b)
128 float_t m = std::max(a,b);
129 return m + std::log(std::exp(a-m) + std::exp(b-m));
145 template<
typename float_t>
146 float_t evalUnivNorm(float_t x, float_t mu, float_t sigma,
bool log)
148 float_t exponent = -.5*(x - mu)*(x-mu)/(sigma*sigma);
151 return -std::log(sigma) - .5*log_two_pi<float_t> + exponent;
153 return inv_sqrt_2pi<float_t> * std::exp(exponent) / sigma;
157 return -std::numeric_limits<float_t>::infinity();
173 template<
typename float_t>
174 float_t evalUnivNorm_unnorm(float_t x, float_t mu, float_t sigma,
bool log)
176 float_t exponent = -.5*(x - mu)*(x-mu)/(sigma*sigma);
181 return std::exp(exponent);
185 return -std::numeric_limits<float_t>::infinity();
198 template<
typename float_t>
199 float_t evalUnivStdNormCDF(float_t x)
202 float_t a1 = 0.254829592;
203 float_t a2 = -0.284496736;
204 float_t a3 = 1.421413741;
205 float_t a4 = -1.453152027;
206 float_t a5 = 1.061405429;
207 float_t p = 0.3275911;
213 float_t xt = std::fabs(x)/std::sqrt(2.0);
216 float_t t = 1.0/(1.0 + p*xt);
217 float_t y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*std::exp(-xt*xt);
219 return 0.5*(1.0 + sign*y);
231 template<
typename float_t>
232 float_t evalUnivBeta(float_t x, float_t alpha, float_t beta,
bool log)
234 if( (x > 0.0) && (x < 1.0) && (alpha > 0.0) && (beta > 0.0) ){
236 return std::lgamma(alpha + beta) - std::lgamma(alpha) - std::lgamma(beta) + (alpha - 1.0)*std::log(x) + (beta - 1.0) * std::log(1.0 - x);
238 return pow(x, alpha-1.0) * pow(1.0-x, beta-1.0) * std::tgamma(alpha + beta) / ( std::tgamma(alpha) * std::tgamma(beta) );
243 return -std::numeric_limits<float_t>::infinity();
259 template<
typename float_t>
260 float_t evalUnivBeta_unnorm(float_t x, float_t alpha, float_t beta,
bool log)
262 if( (x > 0.0) && (x < 1.0) && (alpha > 0.0) && (beta > 0.0) ){
264 return (alpha - 1.0)*std::log(x) + (beta - 1.0) * std::log(1.0 - x);
266 return pow(x, alpha-1.0) * pow(1.0-x, beta-1.0);
271 return -std::numeric_limits<float_t>::infinity();
287 template<
typename float_t>
288 float_t evalUnivInvGamma(float_t x, float_t alpha, float_t beta,
bool log)
290 if ( (x > 0.0) && (alpha > 0.0) && (beta > 0.0) ){
292 return alpha * std::log(beta) - std::lgamma(alpha) - (alpha + 1.0)*std::log(x) - beta/x;
294 return pow(x, -alpha-1.0) * exp(-beta/x) * pow(beta, alpha) / std::tgamma(alpha);
298 return -std::numeric_limits<float_t>::infinity();
314 template<
typename float_t>
315 float_t evalUnivInvGamma_unnorm(float_t x, float_t alpha, float_t beta,
bool log)
317 if ( (x > 0.0) && (alpha > 0.0) && (beta > 0.0) ){
319 return (-alpha - 1.0)*std::log(x) - beta/x;
321 return pow(x, -alpha-1.0) * exp(-beta/x);
325 return -std::numeric_limits<float_t>::infinity();
340 template<
typename float_t>
341 float_t evalUnivHalfNorm(float_t x, float_t sigmaSqd,
bool log)
343 if( (x >= 0.0) && (sigmaSqd > 0.0)){
345 return .5*log_two_over_pi<float_t> - .5*std::log(sigmaSqd) - .5*x*x / sigmaSqd;
347 return std::exp(-.5*x*x/sigmaSqd) * sqrt_two_over_pi<float_t> / std::sqrt(sigmaSqd);
351 return -std::numeric_limits<float_t>::infinity();
366 template<
typename float_t>
367 float_t evalUnivHalfNorm_unnorm(float_t x, float_t sigmaSqd,
bool log)
369 if( (x >= 0.0) && (sigmaSqd > 0.0)){
371 return -.5*x*x / sigmaSqd;
373 return std::exp(-.5*x*x/sigmaSqd);
377 return -std::numeric_limits<float_t>::infinity();
395 template<
typename float_t>
396 float_t evalUnivTruncNorm(float_t x, float_t mu, float_t sigma, float_t lower, float_t upper,
bool log)
398 if( (sigma > 0.0) && (lower <= x) & (x <= upper) ){
400 return evalUnivNorm(x, mu, sigma,
true)
401 - std::log( evalUnivStdNormCDF((upper-mu)/sigma) - evalUnivStdNormCDF((lower-mu)/sigma));
403 return evalUnivNorm(x,mu,sigma,
false)
404 / ( evalUnivStdNormCDF((upper-mu)/sigma) - evalUnivStdNormCDF((lower-mu)/sigma) );
409 return -std::numeric_limits<float_t>::infinity();
427 template<
typename float_t>
428 float_t evalUnivTruncNorm_unnorm(float_t x, float_t mu, float_t sigma, float_t lower, float_t upper,
bool log)
430 if( (sigma > 0.0) && (lower <= x) & (x <= upper) ){
432 return evalUnivNorm_unnorm(x, mu, sigma,
true);
434 return evalUnivNorm_unnorm(x, mu, sigma,
false);
438 return -std::numeric_limits<float_t>::infinity();
454 template<
typename float_t>
455 float_t evalLogitNormal(float_t x, float_t mu, float_t sigma,
bool log)
457 if( (x >= 0.0) && (x <= 1.0) && (sigma > 0.0)){
459 float_t exponent = -.5*(logit(x) - mu)*(logit(x) - mu) / (sigma*sigma);
461 return -std::log(sigma) - .5*log_two_pi<float_t> - std::log(x) - std::log(1.0-x) + exponent;
463 return inv_sqrt_2pi<float_t> * std::exp(exponent) / (x * (1.0-x) * sigma);
467 return -std::numeric_limits<float_t>::infinity();
483 template<
typename float_t>
484 float_t evalLogitNormal_unnorm(float_t x, float_t mu, float_t sigma,
bool log)
486 if( (x >= 0.0) && (x <= 1.0) && (sigma > 0.0)){
488 float_t exponent = -.5*(logit(x) - mu)*(logit(x) - mu) / (sigma*sigma);
490 return -std::log(x) - std::log(1.0-x) + exponent;
492 return std::exp(exponent) / x / (1.0-x);
496 return -std::numeric_limits<float_t>::infinity();
512 template<
typename float_t>
513 float_t evalTwiceFisherNormal(float_t x, float_t mu, float_t sigma,
bool log)
517 if( (x >= -1.0) && (x <= 1.0) && (sigma > 0.0)){
519 float_t exponent = std::log((1.0+x)/(1.0-x)) - mu;
520 exponent = -.5*exponent*exponent/sigma/sigma;
522 return -std::log(sigma) - .5*log_two_pi<float_t> + std::log(2.0) - std::log(1.0+x) - std::log(1.0-x) + exponent;
524 return inv_sqrt_2pi<float_t> * 2.0 * std::exp(exponent)/( (1.0-x)*(1.0+x)*sigma );
528 return -std::numeric_limits<float_t>::infinity();
544 template<
typename float_t>
545 float_t evalTwiceFisherNormal_unnorm(float_t x, float_t mu, float_t sigma,
bool log)
549 if( (x >= -1.0) && (x <= 1.0) && (sigma > 0.0)){
551 float_t exponent = std::log((1.0+x)/(1.0-x)) - mu;
552 exponent = -.5*exponent*exponent/sigma/sigma;
554 return -std::log(1.0+x) - std::log(1.0-x) + exponent;
556 return std::exp(exponent)/(1.0-x)/(1.0+x);
560 return -std::numeric_limits<float_t>::infinity();
576 template<
typename float_t>
577 float_t evalLogNormal(float_t x, float_t mu, float_t sigma,
bool log)
579 if( (x > 0.0) && (sigma > 0.0)){
581 float_t exponent = std::log(x)-mu;
582 exponent = -.5 * exponent * exponent / sigma / sigma;
584 return -std::log(x) - std::log(sigma) - .5*log_two_pi<float_t> + exponent;
586 return inv_sqrt_2pi<float_t> * std::exp(exponent)/(sigma*x);
590 return -std::numeric_limits<float_t>::infinity();
606 template<
typename float_t>
607 float_t evalLogNormal_unnorm(float_t x, float_t mu, float_t sigma,
bool log)
609 if( (x > 0.0) && (sigma > 0.0)){
611 float_t exponent = std::log(x)-mu;
612 exponent = -.5 * exponent * exponent / sigma / sigma;
614 return -std::log(x) + exponent;
616 return std::exp(exponent)/x;
620 return -std::numeric_limits<float_t>::infinity();
636 template<
typename float_t>
637 float_t evalUniform(float_t x, float_t lower, float_t upper,
bool log)
640 if( (x > lower) && (x <= upper)){
642 float_t width = upper-lower;
644 return -std::log(width);
650 return -std::numeric_limits<float_t>::infinity();
667 template<
typename float_t>
668 float_t evalUniform_unnorm(float_t x, float_t lower, float_t upper,
bool log)
671 if( (x > lower) && (x <= upper)){
680 return -std::numeric_limits<float_t>::infinity();
698 template<
typename float_t>
699 float_t evalScaledT(float_t x, float_t mu, float_t sigma, float_t dof,
bool log)
702 if( (sigma > 0.0) && (dof > 0.0) ){
704 float_t zscore = (x-mu)/sigma;
705 float_t lmt = - .5*(dof+1.0)*std::log(1.0 + (zscore*zscore)/dof);
707 return std::lgamma(.5*(dof+1.0)) - std::log(sigma) - .5*std::log(dof) - .5*log_pi<float_t> - std::lgamma(.5*dof) + lmt;
709 return std::exp(std::lgamma(.5*(dof+1.0)) - std::log(sigma) - .5*std::log(dof) - .5*log_pi<float_t> - std::lgamma(.5*dof) + lmt);
712 return -std::numeric_limits<float_t>::infinity();
728 template<
typename float_t>
729 float_t evalScaledT_unnorm(float_t x, float_t mu, float_t sigma, float_t dof,
bool log)
731 if( (sigma > 0.0) && (dof > 0.0) ){
733 float_t zscore = (x-mu)/sigma;
734 float_t lmt = - .5*(dof+1.0)*std::log(1.0 + (zscore*zscore)/dof);
738 return std::exp(lmt);
741 return -std::numeric_limits<float_t>::infinity();
755 template<
typename int_t,
typename float_t>
756 float_t evalDiscreteUnif(int_t x,
int k,
bool log)
758 if( (1 <= x) && (x <= k) ){
760 return -std::log(static_cast<float_t>(k));
762 return 1.0 /
static_cast<float_t
>(k);
766 return -std::numeric_limits<float_t>::infinity();
781 template<
typename int_t,
typename float_t>
782 float_t evalDiscreteUnif_unnorm(int_t x, int_t k,
bool log)
784 if( (1 <= x) && (x <= k) ){
792 return -std::numeric_limits<float_t>::infinity();
806 template<
typename int_t,
typename float_t>
807 float_t evalBernoulli(int_t x, float_t p,
bool log)
809 if( ((x == 0) || (x == 1)) && ( (0.0 <= p) && (p <= 1.0) ) ){
811 return (x==1) ? std::log(p) : std::log(1.0-p);
813 return (x==1) ? p : (1.0-p);
817 return -std::numeric_limits<float_t>::infinity();
833 template<
typename int_t,
typename float_t>
834 float_t evalSkellam(int_t x, float_t mu1, float_t mu2,
bool log)
836 if( (mu1 > 0) && (mu2 > 0) ){
842 using boost::math::tools::evaluate_polynomial;
843 float_t z = 2*std::sqrt(mu1*mu2);
844 float_t log_I(-std::numeric_limits<float_t>::infinity());
863 static const float_t P[]
864 = {1.00000000000000000e+00, 2.49999999999999909e-01,
865 2.77777777777782257e-02, 1.73611111111023792e-03,
866 6.94444444453352521e-05, 1.92901234513219920e-06,
867 3.93675991102510739e-08, 6.15118672704439289e-10,
868 7.59407002058973446e-12, 7.59389793369836367e-14,
869 6.27767773636292611e-16, 4.34709704153272287e-18,
870 2.63417742690109154e-20, 1.13943037744822825e-22,
871 9.07926920085624812e-25};
872 log_I = log_sum_exp<float_t>(0.0, 2.0*std::log(z) - std::log(4.0) + std::log(evaluate_polynomial(P, mu1*mu2)));
878 static const float_t P[]
879 = {3.98942280401425088e-01, 4.98677850604961985e-02,
880 2.80506233928312623e-02, 2.92211225166047873e-02,
881 4.44207299493659561e-02, 1.30970574605856719e-01,
882 -3.35052280231727022e+00, 2.33025711583514727e+02,
883 -1.13366350697172355e+04, 4.24057674317867331e+05,
884 -1.23157028595698731e+07, 2.80231938155267516e+08,
885 -5.01883999713777929e+09, 7.08029243015109113e+10,
886 -7.84261082124811106e+11, 6.76825737854096565e+12,
887 -4.49034849696138065e+13, 2.24155239966958995e+14,
888 -8.13426467865659318e+14, 2.02391097391687777e+15,
889 -3.08675715295370878e+15, 2.17587543863819074e+15};
891 log_I = z + std::log(evaluate_polynomial(P, 1.0/z)) - 0.5*std::log(z);
897 static const float_t P[] = {3.98942280401432905e-01, 4.98677850491434560e-02,
898 2.80506308916506102e-02, 2.92179096853915176e-02,
899 4.53371208762579442e-02};
900 log_I = z + std::log(evaluate_polynomial(P, 1.0/z)) - 0.5*std::log(z);
904 }
else if(std::abs(x)==1){
913 static const float_t P[]
914 = {8.333333333333333803e-02, 6.944444444444341983e-03,
915 3.472222222225921045e-04, 1.157407407354987232e-05,
916 2.755731926254790268e-07, 4.920949692800671435e-09,
917 6.834657311305621830e-11, 7.593969849687574339e-13,
918 6.904822652741917551e-15, 5.220157095351373194e-17,
919 3.410720494727771276e-19, 1.625212890947171108e-21,
920 1.332898928162290861e-23};
923 float_t Q[3] = {1, 0.5, evaluate_polynomial(P, a)};
924 log_I = std::log(z) + std::log(evaluate_polynomial(Q, a)) - std::log(2.0);
930 static const double P[]
931 = {3.989422804014406054e-01, -1.496033551613111533e-01,
932 -4.675104253598537322e-02, -4.090895951581637791e-02,
933 -5.719036414430205390e-02, -1.528189554374492735e-01,
934 3.458284470977172076e+00, -2.426181371595021021e+02,
935 1.178785865993440669e+04, -4.404655582443487334e+05,
936 1.277677779341446497e+07, -2.903390398236656519e+08,
937 5.192386898222206474e+09, -7.313784438967834057e+10,
938 8.087824484994859552e+11, -6.967602516005787001e+12,
939 4.614040809616582764e+13, -2.298849639457172489e+14,
940 8.325554073334618015e+14, -2.067285045778906105e+15,
941 3.146401654361325073e+15, -2.213318202179221945e+15};
942 log_I = z + std::log(evaluate_polynomial(P, 1.0/z)) - 0.5*std::log(z);
947 static const double P[]
948 = {3.989422804014314820e-01, -1.496033551467584157e-01,
949 -4.675105322571775911e-02, -4.090421597376992892e-02,
950 -5.843630344778927582e-02};
951 log_I = z + std::log(evaluate_polynomial(P, 1.0/z)) - 0.5*std::log(z);
957 float_t lim = std::pow((x*x + 2.5) / (2 * z), 3) / 24;
959 if (lim < std::numeric_limits<float_t>::epsilon()* 10) {
963 float_t num = mu - 1;
972 log_I = z - .5*std::log(z) - .5*log_two_pi<float_t> + std::log(s);
978 int_t absx = (x < 0) ? -x : x;
979 float_t lm1m2 = std::log(mu1) + std::log(mu2);
980 float_t first = .5*absx*lm1m2;
981 float_t second = 0.0;
982 float_t third = std::lgamma(absx+1);
983 log_I = first - second - third;
987 float_t last_iter_log_I;
990 second += std::log(m);
991 third += std::log(m+absx);
992 last_iter_log_I = log_I;
993 log_I = log_sum_exp<float_t>(log_I, first - second - third);
996 std::cout <<
"first, second, third, mu1, mu2: " << first <<
", " << second <<
", " << third <<
", " << mu1 <<
", " << mu2 <<
"\n";
997 }
while(log_I != last_iter_log_I);
1002 float_t log_mass = -mu1 - mu2 + .5*x*(std::log(mu1) - std::log(mu2)) + log_I;
1010 return std::exp(log_mass);
1014 return -std::numeric_limits<float_t>::infinity();
1038 template<std::
size_t dim,
typename float_t>
1039 float_t evalMultivNorm(
const Eigen::Matrix<float_t,dim,1> &x,
1040 const Eigen::Matrix<float_t,dim,1> &meanVec,
1041 const Eigen::Matrix<float_t,dim,dim> &covMat,
1044 using Mat = Eigen::Matrix<float_t,dim,dim>;
1045 Eigen::LLT<Mat> lltM(covMat);
1046 if(lltM.info() == Eigen::NumericalIssue)
return log ? -std::numeric_limits<float_t>::infinity() : 0.0;
1047 Mat L = lltM.matrixL();
1048 float_t quadform = (lltM.solve(x-meanVec)).squaredNorm();
1051 for(
size_t i = 0; i < dim; ++i){
1052 ld += std::log(L(i,i));
1056 float_t logDens = -.5*log_two_pi<float_t> * dim - .5*ld - .5*quadform;
1061 return std::exp(logDens);
1078 template<std::
size_t dim,
typename float_t>
1079 float_t evalMultivT(
const Eigen::Matrix<float_t,dim,1> &x,
1080 const Eigen::Matrix<float_t,dim,1> &locVec,
1081 const Eigen::Matrix<float_t,dim,dim> &shapeMat,
1085 if(dof <= 0.0)
return log ? -std::numeric_limits<float_t>::infinity() : 0.0;
1086 using Mat = Eigen::Matrix<float_t,dim,dim>;
1087 Eigen::LLT<Mat> lltM(shapeMat);
1088 if(lltM.info() == Eigen::NumericalIssue)
return log ? -std::numeric_limits<float_t>::infinity() : 0.0;
1089 Mat L = lltM.matrixL();
1090 float_t quadform = (lltM.solve(x-locVec)).squaredNorm();
1093 for(
size_t i = 0; i < dim; ++i){
1094 ld += std::log(L(i,i));
1098 float_t logDens = std::lgamma(.5*(dof+dim)) - .5*dim*std::log(dof) - .5*dim*log_pi<float_t>
1099 - std::lgamma(.5*dof) -.5*ld - .5*(dof+dim)*std::log( 1.0 + quadform/dof );
1104 return std::exp(logDens);
1120 template<std::
size_t bigd, std::
size_t smalld,
typename float_t>
1121 float_t evalMultivNormWBDA(
const Eigen::Matrix<float_t,bigd,1> &x,
1122 const Eigen::Matrix<float_t,bigd,1> &meanVec,
1123 const Eigen::Matrix<float_t,bigd,1> &A,
1124 const Eigen::Matrix<float_t,bigd,smalld> &U,
1125 const Eigen::Matrix<float_t,smalld,smalld> &C,
1129 using bigmat = Eigen::Matrix<float_t,bigd,bigd>;
1130 using smallmat = Eigen::Matrix<float_t,smalld,smalld>;
1132 bigmat Ainv = A.asDiagonal().inverse();
1133 smallmat Cinv = C.inverse();
1134 smallmat I = Cinv + U.transpose()*Ainv*U;
1135 bigmat SigInv = Ainv - Ainv * U * I.ldlt().solve(U.transpose() * Ainv);
1136 Eigen::LLT<bigmat> lltSigInv(SigInv);
1137 bigmat L = lltSigInv.matrixL();
1138 float_t quadform = (L * (x-meanVec)).squaredNorm();
1142 float_t halfld (0.0);
1144 for(
size_t i = 0; i < bigd; ++i){
1145 halfld += std::log(L(i,i));
1148 return -.5*log_two_pi<float_t> * bigd + halfld - .5*quadform;
1152 float_t normConst = std::pow(inv_sqrt_2pi<float_t>, bigd) * L.determinant();
1153 return normConst * std::exp(-.5* quadform);
1171 template<std::
size_t dim,
typename float_t>
1172 float_t evalWishart(
const Eigen::Matrix<float_t,dim,dim> &X,
1173 const Eigen::Matrix<float_t,dim,dim> &Vinv,
1174 const unsigned int &n,
1177 using Mat = Eigen::Matrix<float_t,dim,dim>;
1178 Eigen::LLT<Mat> lltX(X);
1179 Eigen::LLT<Mat> lltVinv(Vinv);
1180 if((n < dim) | (lltX.info() == Eigen::NumericalIssue) | (lltVinv.info() == Eigen::NumericalIssue))
return log ? -std::numeric_limits<float_t>::infinity() : 0.0;
1184 float_t ldvinv (0.0);
1185 Mat Lx = lltX.matrixL();
1186 Mat Lvi = lltVinv.matrixL();
1187 float_t logGammaNOver2 = .25*dim*(dim-1)*log_pi<float_t>;
1190 for(
size_t i = 0; i < dim; ++i){
1191 ldx += std::log(Lx(i,i));
1192 ldvinv += std::log(Lvi(i,i));
1193 logGammaNOver2 += std::lgamma(.5*(n-i));
1198 float_t logDens = .5*(n - dim -1)*ldx - .5*(Vinv*X).trace() - .5*n*dim*std::log(2.0) + .5*n*ldvinv - logGammaNOver2;
1203 return std::exp(logDens);
1220 template<std::
size_t dim,
typename float_t>
1221 float_t evalInvWishart(
const Eigen::Matrix<float_t,dim,dim> &X,
1222 const Eigen::Matrix<float_t,dim,dim> &Psi,
1223 const unsigned int &nu,
1226 using Mat = Eigen::Matrix<float_t,dim,dim>;
1227 Eigen::LLT<Mat> lltX(X);
1228 Eigen::LLT<Mat> lltPsi(Psi);
1229 if((nu < dim) | (lltX.info() == Eigen::NumericalIssue) | (lltPsi.info() == Eigen::NumericalIssue))
return log ? -std::numeric_limits<float_t>::infinity() : 0.0;
1233 float_t ldPsi (0.0);
1234 Mat Lx = lltX.matrixL();
1235 Mat Lpsi = lltPsi.matrixL();
1236 float_t logGammaNuOver2 = .25*dim*(dim-1)*log_pi<float_t>;
1239 for(
size_t i = 0; i < dim; ++i){
1240 ldx += std::log(Lx(i,i));
1241 ldPsi += std::log(Lpsi(i,i));
1242 logGammaNuOver2 += std::lgamma(.5*(nu-i));
1249 float_t logDens = .5*nu*ldPsi - .5*(nu + dim + 1.0)*ldx - .5*(Psi*X.inverse() ).trace() - .5*nu*dim*std::log(2.0) - logGammaNuOver2;
1255 return std::exp(logDens);