pf
rv_eval.h
1 #ifndef RV_EVAL_H
2 #define RV_EVAL_H
3 
4 #include <cstddef> // std::size_t
5 
6 #ifdef DROPPINGTHISINRPACKAGE
7  #include <RcppEigen.h>
8  // [[Rcpp::depends(RcppEigen)]]
9 #else
10  #include <Eigen/Dense>
11 #endif
12 
13 #include <iostream> // cerr
14 #include "boost/math/special_functions.hpp"
15 
16 
17 namespace pf {
18 
19 namespace rveval{
20 
24 
26 template<class T>
27 constexpr T inv_sqrt_2pi = T(0.3989422804014327);
28 
30 template<class T>
31 constexpr T sqrt_two_over_pi(0.797884560802865);
32 
34 template<class T>
35 constexpr T log_two_pi (1.83787706640935);
36 
38 template<class T>
39 constexpr T log_two_over_pi (-0.451582705289455);
40 
42 template<class T>
43 constexpr T log_pi (1.1447298858494);
44 
45 
49 
55 template<typename float_t>
56 float_t twiceFisher(float_t phi)
57 {
58  if ( (phi <= -1.0) || (phi >= 1.0) )
59  throw std::invalid_argument( "error: phi was not between -1 and 1" );
60  else
61  return std::log(1.0 + phi) - std::log(1.0 - phi);
62 }
63 
64 
65 
71 template<typename float_t>
72 float_t invTwiceFisher(float_t psi)
73 {
74  float_t ans = (1.0 - std::exp(psi)) / ( -1.0 - std::exp(psi) );
75 
76  if ( (ans <= -1.0) || (ans >= 1.0) )
77  std::cerr << "error: there was probably overflow for exp(psi) \n";
78 
79  return ans;
80 }
81 
82 
88 template<typename float_t>
89 float_t logit(float_t p)
90 {
91  if ( (p <= 0.0) || (p >= 1.0))
92  std::cerr << "error: p was not between 0 and 1 \n";
93 
94  return std::log(p) - std::log(1.0 - p);
95 }
96 
97 
103 template<typename float_t>
104 float_t inv_logit(float_t r)
105 {
106  float_t ans = 1.0/( 1.0 + std::exp(-r) );
107 
108  if ( (ans <= 0.0) || (ans >= 1.0))
109  std::cerr << "error: there was probably underflow for exp(-r) \n";
110 
111  return ans;
112 }
113 
114 
120 template<typename float_t>
121 float_t log_inv_logit(float_t r)
122 {
123  if(r < -750.00 || r > 750.00) std::cerr << "warning: log_inv_logit might be under/over-flowing\n";
124  return -std::log(1.0 + std::exp(-r));
125 }
126 
127 
134 template<typename float_t>
135 float_t log_sum_exp(float_t a, float_t b)
136 {
137  float_t m = std::max(a,b);
138  return m + std::log(std::exp(a-m) + std::exp(b-m));
139 }
140 
141 
145 
154 template<typename float_t>
155 float_t evalUnivNorm(float_t x, float_t mu, float_t sigma, bool log)
156 {
157  float_t exponent = -.5*(x - mu)*(x-mu)/(sigma*sigma);
158  if( sigma > 0.0){
159  if(log){
160  return -std::log(sigma) - .5*log_two_pi<float_t> + exponent;
161  }else{
162  return inv_sqrt_2pi<float_t> * std::exp(exponent) / sigma;
163  }
164  }else{
165  if(log){
166  return -std::numeric_limits<float_t>::infinity();
167  }else{
168  return 0.0;
169  }
170  }
171 }
172 
173 
182 template<typename float_t>
183 float_t evalUnivNorm_unnorm(float_t x, float_t mu, float_t sigma, bool log)
184 {
185  float_t exponent = -.5*(x - mu)*(x-mu)/(sigma*sigma);
186  if( sigma > 0.0){
187  if(log){
188  return exponent;
189  }else{
190  return std::exp(exponent);
191  }
192  }else{
193  if(log){
194  return -std::numeric_limits<float_t>::infinity();
195  }else{
196  return 0.0;
197  }
198  }
199 }
200 
201 
207 template<typename float_t>
208 float_t evalUnivStdNormCDF(float_t x) // john cook code
209 {
210  // constants
211  float_t a1 = 0.254829592;
212  float_t a2 = -0.284496736;
213  float_t a3 = 1.421413741;
214  float_t a4 = -1.453152027;
215  float_t a5 = 1.061405429;
216  float_t p = 0.3275911;
217 
218  // Save the sign of x
219  int sign = 1;
220  if (x < 0)
221  sign = -1;
222  float_t xt = std::fabs(x)/std::sqrt(2.0);
223 
224  // A&S formula 7.1.26
225  float_t t = 1.0/(1.0 + p*xt);
226  float_t y = 1.0 - (((((a5*t + a4)*t) + a3)*t + a2)*t + a1)*t*std::exp(-xt*xt);
227 
228  return 0.5*(1.0 + sign*y);
229 }
230 
231 
240 template<typename float_t>
241 float_t evalUnivBeta(float_t x, float_t alpha, float_t beta, bool log)
242 {
243  if( (x > 0.0) && (x < 1.0) && (alpha > 0.0) && (beta > 0.0) ){ // x in support and parameters acceptable
244  if(log){
245  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);
246  }else{
247  return pow(x, alpha-1.0) * pow(1.0-x, beta-1.0) * std::tgamma(alpha + beta) / ( std::tgamma(alpha) * std::tgamma(beta) );
248  }
249 
250  }else{ //not ( x in support and parameters acceptable )
251  if(log){
252  return -std::numeric_limits<float_t>::infinity();
253  }else{
254  return 0.0;
255  }
256  }
257 }
258 
259 
268 template<typename float_t>
269 float_t evalUnivBeta_unnorm(float_t x, float_t alpha, float_t beta, bool log)
270 {
271  if( (x > 0.0) && (x < 1.0) && (alpha > 0.0) && (beta > 0.0) ){ // x in support and parameters acceptable
272  if(log){
273  return (alpha - 1.0)*std::log(x) + (beta - 1.0) * std::log(1.0 - x);
274  }else{
275  return pow(x, alpha-1.0) * pow(1.0-x, beta-1.0);
276  }
277 
278  }else{ //not ( x in support and parameters acceptable )
279  if(log){
280  return -std::numeric_limits<float_t>::infinity();
281  }else{
282  return 0.0;
283  }
284  }
285 }
286 
287 
296 template<typename float_t>
297 float_t evalUnivInvGamma(float_t x, float_t alpha, float_t beta, bool log)
298 {
299  if ( (x > 0.0) && (alpha > 0.0) && (beta > 0.0) ){ // x in support and acceptable parameters
300  if (log){
301  return alpha * std::log(beta) - std::lgamma(alpha) - (alpha + 1.0)*std::log(x) - beta/x;
302  }else{
303  return pow(x, -alpha-1.0) * exp(-beta/x) * pow(beta, alpha) / std::tgamma(alpha);
304  }
305  }else{ // not ( x in support and acceptable parameters )
306  if (log){
307  return -std::numeric_limits<float_t>::infinity();
308  }else{
309  return 0.0;
310  }
311  }
312 }
313 
314 
323 template<typename float_t>
324 float_t evalUnivInvGamma_unnorm(float_t x, float_t alpha, float_t beta, bool log)
325 {
326  if ( (x > 0.0) && (alpha > 0.0) && (beta > 0.0) ){ // x in support and acceptable parameters
327  if (log){
328  return (-alpha - 1.0)*std::log(x) - beta/x;
329  }else{
330  return pow(x, -alpha-1.0) * exp(-beta/x);
331  }
332  }else{ // not ( x in support and acceptable parameters )
333  if (log){
334  return -std::numeric_limits<float_t>::infinity();
335  }else{
336  return 0.0;
337  }
338  }
339 }
340 
341 
349 template<typename float_t>
350 float_t evalUnivHalfNorm(float_t x, float_t sigmaSqd, bool log)
351 {
352  if( (x >= 0.0) && (sigmaSqd > 0.0)){
353  if (log){
354  return .5*log_two_over_pi<float_t> - .5*std::log(sigmaSqd) - .5*x*x / sigmaSqd;
355  }else{
356  return std::exp(-.5*x*x/sigmaSqd) * sqrt_two_over_pi<float_t> / std::sqrt(sigmaSqd);
357  }
358  }else{
359  if (log){
360  return -std::numeric_limits<float_t>::infinity();
361  }else{
362  return 0.0;
363  }
364  }
365 }
366 
367 
375 template<typename float_t>
376 float_t evalUnivHalfNorm_unnorm(float_t x, float_t sigmaSqd, bool log)
377 {
378  if( (x >= 0.0) && (sigmaSqd > 0.0)){
379  if (log){
380  return -.5*x*x / sigmaSqd;
381  }else{
382  return std::exp(-.5*x*x/sigmaSqd);
383  }
384  }else{
385  if (log){
386  return -std::numeric_limits<float_t>::infinity();
387  }else{
388  return 0.0;
389  }
390  }
391 }
392 
393 
404 template<typename float_t>
405 float_t evalUnivTruncNorm(float_t x, float_t mu, float_t sigma, float_t lower, float_t upper, bool log)
406 {
407  if( (sigma > 0.0) && (lower <= x) & (x <= upper) ){
408  if(log){
409  return evalUnivNorm(x, mu, sigma, true)
410  - std::log( evalUnivStdNormCDF((upper-mu)/sigma) - evalUnivStdNormCDF((lower-mu)/sigma));
411  }else{
412  return evalUnivNorm(x,mu,sigma,false)
413  / ( evalUnivStdNormCDF((upper-mu)/sigma) - evalUnivStdNormCDF((lower-mu)/sigma) );
414  }
415 
416  }else{
417  if (log){
418  return -std::numeric_limits<float_t>::infinity();
419  }else{
420  return 0.0;
421  }
422  }
423 }
424 
425 
436 template<typename float_t>
437 float_t evalUnivTruncNorm_unnorm(float_t x, float_t mu, float_t sigma, float_t lower, float_t upper, bool log)
438 {
439  if( (sigma > 0.0) && (lower <= x) & (x <= upper) ){
440  if(log){
441  return evalUnivNorm_unnorm(x, mu, sigma, true);
442  }else{
443  return evalUnivNorm_unnorm(x, mu, sigma, false);
444  }
445  }else{
446  if (log){
447  return -std::numeric_limits<float_t>::infinity();
448  }else{
449  return 0.0;
450  }
451  }
452 }
453 
454 
463 template<typename float_t>
464 float_t evalLogitNormal(float_t x, float_t mu, float_t sigma, bool log)
465 {
466  if( (x >= 0.0) && (x <= 1.0) && (sigma > 0.0)){
467 
468  float_t exponent = -.5*(logit(x) - mu)*(logit(x) - mu) / (sigma*sigma);
469  if(log){
470  return -std::log(sigma) - .5*log_two_pi<float_t> - std::log(x) - std::log(1.0-x) + exponent;
471  }else{
472  return inv_sqrt_2pi<float_t> * std::exp(exponent) / (x * (1.0-x) * sigma);
473  }
474  }else{
475  if(log){
476  return -std::numeric_limits<float_t>::infinity();
477  }else{
478  return 0.0;
479  }
480  }
481 }
482 
483 
492 template<typename float_t>
493 float_t evalLogitNormal_unnorm(float_t x, float_t mu, float_t sigma, bool log)
494 {
495  if( (x >= 0.0) && (x <= 1.0) && (sigma > 0.0)){
496 
497  float_t exponent = -.5*(logit(x) - mu)*(logit(x) - mu) / (sigma*sigma);
498  if(log){
499  return -std::log(x) - std::log(1.0-x) + exponent;
500  }else{
501  return std::exp(exponent) / x / (1.0-x);
502  }
503  }else{
504  if(log){
505  return -std::numeric_limits<float_t>::infinity();
506  }else{
507  return 0.0;
508  }
509  }
510 }
511 
512 
521 template<typename float_t>
522 float_t evalTwiceFisherNormal(float_t x, float_t mu, float_t sigma, bool log)
523 {
524 
525  // https://stats.stackexchange.com/questions/321905/what-is-the-name-of-this-random-variable/321907#321907
526  if( (x >= -1.0) && (x <= 1.0) && (sigma > 0.0)){
527 
528  float_t exponent = std::log((1.0+x)/(1.0-x)) - mu;
529  exponent = -.5*exponent*exponent/sigma/sigma;
530  if(log){
531  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;
532  }else{
533  return inv_sqrt_2pi<float_t> * 2.0 * std::exp(exponent)/( (1.0-x)*(1.0+x)*sigma );
534  }
535  }else{
536  if(log){
537  return -std::numeric_limits<float_t>::infinity();
538  }else{
539  return 0.0;
540  }
541  }
542 }
543 
544 
553 template<typename float_t>
554 float_t evalTwiceFisherNormal_unnorm(float_t x, float_t mu, float_t sigma, bool log)
555 {
556 
557  // https://stats.stackexchange.com/questions/321905/what-is-the-name-of-this-random-variable/321907#321907
558  if( (x >= -1.0) && (x <= 1.0) && (sigma > 0.0)){
559 
560  float_t exponent = std::log((1.0+x)/(1.0-x)) - mu;
561  exponent = -.5*exponent*exponent/sigma/sigma;
562  if(log){
563  return -std::log(1.0+x) - std::log(1.0-x) + exponent;
564  }else{
565  return std::exp(exponent)/(1.0-x)/(1.0+x);
566  }
567  }else{
568  if(log){
569  return -std::numeric_limits<float_t>::infinity();
570  }else{
571  return 0.0;
572  }
573  }
574 }
575 
576 
585 template<typename float_t>
586 float_t evalLogNormal(float_t x, float_t mu, float_t sigma, bool log)
587 {
588  if( (x > 0.0) && (sigma > 0.0)){
589 
590  float_t exponent = std::log(x)-mu;
591  exponent = -.5 * exponent * exponent / sigma / sigma;
592  if(log){
593  return -std::log(x) - std::log(sigma) - .5*log_two_pi<float_t> + exponent;
594  }else{
595  return inv_sqrt_2pi<float_t> * std::exp(exponent)/(sigma*x);
596  }
597  }else{
598  if(log){
599  return -std::numeric_limits<float_t>::infinity();
600  }else{
601  return 0.0;
602  }
603  }
604 }
605 
606 
615 template<typename float_t>
616 float_t evalLogNormal_unnorm(float_t x, float_t mu, float_t sigma, bool log)
617 {
618  if( (x > 0.0) && (sigma > 0.0)){
619 
620  float_t exponent = std::log(x)-mu;
621  exponent = -.5 * exponent * exponent / sigma / sigma;
622  if(log){
623  return -std::log(x) + exponent;
624  }else{
625  return std::exp(exponent)/x;
626  }
627  }else{
628  if(log){
629  return -std::numeric_limits<float_t>::infinity();
630  }else{
631  return 0.0;
632  }
633  }
634 }
635 
636 
645 template<typename float_t>
646 float_t evalUniform(float_t x, float_t lower, float_t upper, bool log)
647 {
648 
649  if( (x > lower) && (x <= upper)){
650 
651  float_t width = upper-lower;
652  if(log){
653  return -std::log(width);
654  }else{
655  return 1.0/width;
656  }
657  }else{
658  if(log){
659  return -std::numeric_limits<float_t>::infinity();
660  }else{
661  return 0.0;
662  }
663  }
664 
665 }
666 
667 
676 template<typename float_t>
677 float_t evalUniform_unnorm(float_t x, float_t lower, float_t upper, bool log)
678 {
679 
680  if( (x > lower) && (x <= upper)){
681 
682  if(log){
683  return 0.0;
684  }else{
685  return 1.0;
686  }
687  }else{
688  if(log){
689  return -std::numeric_limits<float_t>::infinity();
690  }else{
691  return 0.0;
692  }
693  }
694 
695 }
696 
697 
707 template<typename float_t>
708 float_t evalScaledT(float_t x, float_t mu, float_t sigma, float_t dof, bool log)
709 {
710 
711  if( (sigma > 0.0) && (dof > 0.0) ){
712 
713  float_t zscore = (x-mu)/sigma;
714  float_t lmt = - .5*(dof+1.0)*std::log(1.0 + (zscore*zscore)/dof);
715  if(log)
716  return std::lgamma(.5*(dof+1.0)) - std::log(sigma) - .5*std::log(dof) - .5*log_pi<float_t> - std::lgamma(.5*dof) + lmt;
717  else
718  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);
719  } else{
720  if(log)
721  return -std::numeric_limits<float_t>::infinity();
722  else
723  return 0.0;
724  }
725 }
726 
727 
737 template<typename float_t>
738 float_t evalScaledT_unnorm(float_t x, float_t mu, float_t sigma, float_t dof, bool log)
739 {
740  if( (sigma > 0.0) && (dof > 0.0) ){
741 
742  float_t zscore = (x-mu)/sigma;
743  float_t lmt = - .5*(dof+1.0)*std::log(1.0 + (zscore*zscore)/dof);
744  if(log)
745  return lmt;
746  else
747  return std::exp(lmt);
748  } else{
749  if(log)
750  return -std::numeric_limits<float_t>::infinity();
751  else
752  return 0.0;
753  }
754 }
755 
756 
764 template<typename int_t, typename float_t>
765 float_t evalDiscreteUnif(int_t x, int k, bool log)
766 {
767  if( (1 <= x) && (x <= k) ){
768  if(log){
769  return -std::log(static_cast<float_t>(k));
770  }else{
771  return 1.0 / static_cast<float_t>(k);
772  }
773  }else{ // x not in support
774  if(log){
775  return -std::numeric_limits<float_t>::infinity();
776  }else{
777  return 0.0;
778  }
779  }
780 }
781 
782 
790 template<typename int_t, typename float_t>
791 float_t evalDiscreteUnif_unnorm(int_t x, int_t k, bool log)
792 {
793  if( (1 <= x) && (x <= k) ){
794  if(log){
795  return 0.0;
796  }else{
797  return 1.0;
798  }
799  }else{
800  if(log){
801  return -std::numeric_limits<float_t>::infinity();
802  }else{
803  return 0.0;
804  }
805  }
806 }
807 
808 
815 template<typename int_t, typename float_t>
816 float_t evalBernoulli(int_t x, float_t p, bool log)
817 {
818  if( ((x == 0) || (x == 1)) && ( (0.0 <= p) && (p <= 1.0) ) ){ // if valid x and valid p
819  if(log){
820  return (x==1) ? std::log(p) : std::log(1.0-p);
821  }else{
822  return (x==1) ? p : (1.0-p);
823  }
824  }else{ // either invalid x or invalid p
825  if(log){
826  return -std::numeric_limits<float_t>::infinity();
827  }else{
828  return 0.0;
829  }
830  }
831 }
832 
833 
842 template<typename int_t, typename float_t>
843 float_t evalSkellam(int_t x, float_t mu1, float_t mu2, bool log)
844 {
845  if( (mu1 > 0) && (mu2 > 0) ){
846 
847  // much of this function is adapted from
848  // https://github.com/stan-dev/math/blob/9b2e93ba58fa00521275b22a190468ab22f744a3/stan/math/prim/fun/log_modified_bessel_first_kind.hpp
849 
850  // step 1: calculate log I_k(2\sqrt{mu_1 mu_2}) using log_sum_exp
851  using boost::math::tools::evaluate_polynomial;
852  float_t z = 2*std::sqrt(mu1*mu2);
853  float_t log_I(-std::numeric_limits<float_t>::infinity());
854 
855  if (x == 0) {
856  // modified from Boost's bessel_i0_imp in the double precision case,
857  // which refers to:
858  // Modified Bessel function of the first kind of order zero
859  // we use the approximating forms derived in:
860  // "Rational Approximations for the Modified Bessel Function of the
861  // First Kind -- I0(x) for Computations with Double Precision"
862  // by Pavel Holoborodko, see
863  // http://www.advanpix.com/2015/11/11/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i0-computations-double-precision
864  // The actual coefficients used are [Boost's] own, and extend
865  // Pavel's work to precisions other than double.
866 
867  if (z < 7.75) {
868  // Bessel I0 over[10 ^ -16, 7.75]
869  // Max error in interpolated form : 3.042e-18
870  // Max Error found at double precision = Poly : 5.106609e-16
871  // Cheb : 5.239199e-16
872  static const float_t P[]
873  = {1.00000000000000000e+00, 2.49999999999999909e-01,
874  2.77777777777782257e-02, 1.73611111111023792e-03,
875  6.94444444453352521e-05, 1.92901234513219920e-06,
876  3.93675991102510739e-08, 6.15118672704439289e-10,
877  7.59407002058973446e-12, 7.59389793369836367e-14,
878  6.27767773636292611e-16, 4.34709704153272287e-18,
879  2.63417742690109154e-20, 1.13943037744822825e-22,
880  9.07926920085624812e-25};
881  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)));
882 
883  }else if (z < 500) {
884  // Max error in interpolated form : 1.685e-16
885  // Max Error found at double precision = Poly : 2.575063e-16
886  // Cheb : 2.247615e+00
887  static const float_t P[]
888  = {3.98942280401425088e-01, 4.98677850604961985e-02,
889  2.80506233928312623e-02, 2.92211225166047873e-02,
890  4.44207299493659561e-02, 1.30970574605856719e-01,
891  -3.35052280231727022e+00, 2.33025711583514727e+02,
892  -1.13366350697172355e+04, 4.24057674317867331e+05,
893  -1.23157028595698731e+07, 2.80231938155267516e+08,
894  -5.01883999713777929e+09, 7.08029243015109113e+10,
895  -7.84261082124811106e+11, 6.76825737854096565e+12,
896  -4.49034849696138065e+13, 2.24155239966958995e+14,
897  -8.13426467865659318e+14, 2.02391097391687777e+15,
898  -3.08675715295370878e+15, 2.17587543863819074e+15};
899 
900  log_I = z + std::log(evaluate_polynomial(P, 1.0/z)) - 0.5*std::log(z);
901 
902  }else{
903 
904  // Max error in interpolated form : 2.437e-18
905  // Max Error found at double precision = Poly : 1.216719e-16
906  static const float_t P[] = {3.98942280401432905e-01, 4.98677850491434560e-02,
907  2.80506308916506102e-02, 2.92179096853915176e-02,
908  4.53371208762579442e-02};
909  log_I = z + std::log(evaluate_polynomial(P, 1.0/z)) - 0.5*std::log(z);
910 
911  }
912 
913  }else if(std::abs(x)==1){
914 
915  // modified from Boost's bessel_i1_imp in the double precision case
916  // see credits above in the v == 0 case
917  if (z < 7.75) {
918  // Bessel I0 over[10 ^ -16, 7.75]
919  // Max error in interpolated form: 5.639e-17
920  // Max Error found at double precision = Poly: 1.795559e-16
921 
922  static const float_t P[]
923  = {8.333333333333333803e-02, 6.944444444444341983e-03,
924  3.472222222225921045e-04, 1.157407407354987232e-05,
925  2.755731926254790268e-07, 4.920949692800671435e-09,
926  6.834657311305621830e-11, 7.593969849687574339e-13,
927  6.904822652741917551e-15, 5.220157095351373194e-17,
928  3.410720494727771276e-19, 1.625212890947171108e-21,
929  1.332898928162290861e-23};
930 
931  float_t a = mu1*mu2;
932  float_t Q[3] = {1, 0.5, evaluate_polynomial(P, a)};
933  log_I = std::log(z) + std::log(evaluate_polynomial(Q, a)) - std::log(2.0);
934 
935  }else if (z < 500) {
936  // Max error in interpolated form: 1.796e-16
937  // Max Error found at double precision = Poly: 2.898731e-16
938 
939  static const double P[]
940  = {3.989422804014406054e-01, -1.496033551613111533e-01,
941  -4.675104253598537322e-02, -4.090895951581637791e-02,
942  -5.719036414430205390e-02, -1.528189554374492735e-01,
943  3.458284470977172076e+00, -2.426181371595021021e+02,
944  1.178785865993440669e+04, -4.404655582443487334e+05,
945  1.277677779341446497e+07, -2.903390398236656519e+08,
946  5.192386898222206474e+09, -7.313784438967834057e+10,
947  8.087824484994859552e+11, -6.967602516005787001e+12,
948  4.614040809616582764e+13, -2.298849639457172489e+14,
949  8.325554073334618015e+14, -2.067285045778906105e+15,
950  3.146401654361325073e+15, -2.213318202179221945e+15};
951  log_I = z + std::log(evaluate_polynomial(P, 1.0/z)) - 0.5*std::log(z);
952  }else {
953 
954  // Max error in interpolated form: 1.320e-19
955  // Max Error found at double precision = Poly: 7.065357e-17
956  static const double P[]
957  = {3.989422804014314820e-01, -1.496033551467584157e-01,
958  -4.675105322571775911e-02, -4.090421597376992892e-02,
959  -5.843630344778927582e-02};
960  log_I = z + std::log(evaluate_polynomial(P, 1.0/z)) - 0.5*std::log(z);
961  }
962 
963  }else if( z > 100){
964 
965  // Boost does something like this in asymptotic_bessel_i_large_x
966  float_t lim = std::pow((x*x + 2.5) / (2 * z), 3) / 24;
967 
968  if (lim < std::numeric_limits<float_t>::epsilon()* 10) {
969  float_t s = 1;
970  float_t mu = 4*x*x;
971  float_t ex = 8 * z;
972  float_t num = mu - 1;
973  float_t denom = ex;
974  s -= num / denom;
975  num *= mu - 9;
976  denom *= ex * 2;
977  s += num / denom;
978  num *= mu - 25;
979  denom *= ex * 3;
980  s -= num / denom;
981  log_I = z - .5*std::log(z) - .5*log_two_pi<float_t> + std::log(s);
982  }
983  }else{
984 
985  // just do the sum straightforwardly
986  // m=0
987  int_t absx = (x < 0) ? -x : x;
988  float_t lm1m2 = std::log(mu1) + std::log(mu2);
989  float_t first = .5*absx*lm1m2;
990  float_t second = 0.0;
991  float_t third = std::lgamma(absx+1);
992  log_I = first - second - third;
993 
994  // m > 0
995  float_t m = 1.0;
996  float_t last_iter_log_I;
997  do{
998  first += lm1m2;
999  second += std::log(m);
1000  third += std::log(m+absx);
1001  last_iter_log_I = log_I;
1002  log_I = log_sum_exp<float_t>(log_I, first - second - third);
1003  m++;
1004  if(m > 1000)
1005  std::cout << "first, second, third, mu1, mu2: " << first << ", " << second << ", " << third << ", " << mu1 << ", " << mu2 << "\n";
1006  }while(log_I != last_iter_log_I);
1007 
1008  }
1009 
1010  // step 2: add the easy parts to get the overall pmf evaluation
1011  float_t log_mass = -mu1 - mu2 + .5*x*(std::log(mu1) - std::log(mu2)) + log_I;
1012  //std::cout << "guaranteed: " << std::log(boost::math::cyl_bessel_i<int_t, float_t>(x,2.0*std::sqrt(mu1*mu2)))
1013  // << "\nmine: " << log_I << "\n";
1014 
1015  // step 3: handle log/nonlog particulars
1016  if(log) {
1017  return log_mass;
1018  }else {
1019  return std::exp(log_mass);
1020  }
1021  }else {
1022  if(log) {
1023  return -std::numeric_limits<float_t>::infinity();
1024  }else {
1025  return 0.0;
1026  }
1027  }
1028 }
1029 
1033 
1034 
1047 template<std::size_t dim, typename float_t>
1048 float_t evalMultivNorm(const Eigen::Matrix<float_t,dim,1> &x,
1049  const Eigen::Matrix<float_t,dim,1> &meanVec,
1050  const Eigen::Matrix<float_t,dim,dim> &covMat,
1051  bool log = false)
1052 {
1053  using Mat = Eigen::Matrix<float_t,dim,dim>;
1054  Eigen::LLT<Mat> lltM(covMat);
1055  if(lltM.info() == Eigen::NumericalIssue) return log ? -std::numeric_limits<float_t>::infinity() : 0.0; // if not pd return 0 dens
1056  Mat L = lltM.matrixL(); // the lower diagonal L such that M = LL^T
1057  float_t quadform = (lltM.solve(x-meanVec)).squaredNorm();
1058  float_t ld (0.0); // calculate log-determinant using cholesky decomposition too
1059  // add up log of diagnols of Cholesky L
1060  for(size_t i = 0; i < dim; ++i){
1061  ld += std::log(L(i,i));
1062  }
1063  ld *= 2; // covMat = LL^T
1064 
1065  float_t logDens = -.5*log_two_pi<float_t> * dim - .5*ld - .5*quadform;
1066 
1067  if(log){
1068  return logDens;
1069  }else{
1070  return std::exp(logDens);
1071  }
1072 }
1073 
1074 
1087 template<std::size_t dim, typename float_t>
1088 float_t evalMultivT(const Eigen::Matrix<float_t,dim,1> &x,
1089  const Eigen::Matrix<float_t,dim,1> &locVec,
1090  const Eigen::Matrix<float_t,dim,dim> &shapeMat,
1091  const float_t& dof,
1092  bool log = false)
1093 {
1094  if(dof <= 0.0) return log ? -std::numeric_limits<float_t>::infinity() : 0.0; // degrees of freedom must be positive
1095  using Mat = Eigen::Matrix<float_t,dim,dim>;
1096  Eigen::LLT<Mat> lltM(shapeMat);
1097  if(lltM.info() == Eigen::NumericalIssue) return log ? -std::numeric_limits<float_t>::infinity() : 0.0; // if not pd return 0 dens
1098  Mat L = lltM.matrixL(); // the lower diagonal L such that M = LL^T
1099  float_t quadform = (lltM.solve(x-locVec)).squaredNorm();
1100  float_t ld (0.0); // calculate log-determinant using cholesky decomposition too
1101  // add up log of diagnols of Cholesky L
1102  for(size_t i = 0; i < dim; ++i){
1103  ld += std::log(L(i,i));
1104  }
1105  ld *= 2; // shapeMat = LL^T
1106 
1107  float_t logDens = std::lgamma(.5*(dof+dim)) - .5*dim*std::log(dof) - .5*dim*log_pi<float_t>
1108  - std::lgamma(.5*dof) -.5*ld - .5*(dof+dim)*std::log( 1.0 + quadform/dof );
1109 
1110  if(log){
1111  return logDens;
1112  }else{
1113  return std::exp(logDens);
1114  }
1115 }
1116 
1117 
1129 template<std::size_t bigd, std::size_t smalld, typename float_t>
1130 float_t evalMultivNormWBDA(const Eigen::Matrix<float_t,bigd,1> &x,
1131  const Eigen::Matrix<float_t,bigd,1> &meanVec,
1132  const Eigen::Matrix<float_t,bigd,1> &A,
1133  const Eigen::Matrix<float_t,bigd,smalld> &U,
1134  const Eigen::Matrix<float_t,smalld,smalld> &C,
1135  bool log = false)
1136 {
1137 
1138  using bigmat = Eigen::Matrix<float_t,bigd,bigd>;
1139  using smallmat = Eigen::Matrix<float_t,smalld,smalld>;
1140 
1141  bigmat Ainv = A.asDiagonal().inverse();
1142  smallmat Cinv = C.inverse();
1143  smallmat I = Cinv + U.transpose()*Ainv*U;
1144  bigmat SigInv = Ainv - Ainv * U * I.ldlt().solve(U.transpose() * Ainv);
1145  Eigen::LLT<bigmat> lltSigInv(SigInv);
1146  bigmat L = lltSigInv.matrixL(); // LL' = Sig^{-1}
1147  float_t quadform = (L * (x-meanVec)).squaredNorm();
1148  if (log){
1149 
1150  // calculate log-determinant using cholesky decomposition (assumes symmetric and positive definite)
1151  float_t halfld (0.0);
1152  // add up log of diagnols of Cholesky L
1153  for(size_t i = 0; i < bigd; ++i){
1154  halfld += std::log(L(i,i));
1155  }
1156 
1157  return -.5*log_two_pi<float_t> * bigd + halfld - .5*quadform;
1158 
1159 
1160  }else{ // not the log density
1161  float_t normConst = std::pow(inv_sqrt_2pi<float_t>, bigd) * L.determinant();
1162  return normConst * std::exp(-.5* quadform);
1163  }
1164 
1165 }
1166 
1167 
1180 template<std::size_t dim, typename float_t>
1181 float_t evalWishart(const Eigen::Matrix<float_t,dim,dim> &X,
1182  const Eigen::Matrix<float_t,dim,dim> &Vinv,
1183  const unsigned int &n,
1184  bool log = false)
1185 {
1186  using Mat = Eigen::Matrix<float_t,dim,dim>;
1187  Eigen::LLT<Mat> lltX(X);
1188  Eigen::LLT<Mat> lltVinv(Vinv);
1189  if((n < dim) | (lltX.info() == Eigen::NumericalIssue) | (lltVinv.info() == Eigen::NumericalIssue)) return log ? -std::numeric_limits<float_t>::infinity() : 0.0;
1190  // https://stackoverflow.com/questions/35227131/eigen-check-if-matrix-is-positive-semi-definite
1191 
1192  float_t ldx (0.0); // log determinant of X
1193  float_t ldvinv (0.0);
1194  Mat Lx = lltX.matrixL(); // the lower diagonal L such that X = LL^T
1195  Mat Lvi = lltVinv.matrixL();
1196  float_t logGammaNOver2 = .25*dim*(dim-1)*log_pi<float_t>; // existence guaranteed when n > dim-1
1197 
1198  // add up log of diagonals of each Cholesky L
1199  for(size_t i = 0; i < dim; ++i){
1200  ldx += std::log(Lx(i,i));
1201  ldvinv += std::log(Lvi(i,i));
1202  logGammaNOver2 += std::lgamma(.5*(n-i)); // recall j = i+1
1203  }
1204  ldx *= 2.0; // X = LL^T
1205  ldvinv *= 2.0;
1206 
1207  float_t logDens = .5*(n - dim -1)*ldx - .5*(Vinv*X).trace() - .5*n*dim*std::log(2.0) + .5*n*ldvinv - logGammaNOver2;
1208 
1209  if(log){
1210  return logDens;
1211  }else{
1212  return std::exp(logDens);
1213  }
1214 }
1215 
1216 
1229 template<std::size_t dim, typename float_t>
1230 float_t evalInvWishart(const Eigen::Matrix<float_t,dim,dim> &X,
1231  const Eigen::Matrix<float_t,dim,dim> &Psi,
1232  const unsigned int &nu,
1233  bool log = false)
1234 {
1235  using Mat = Eigen::Matrix<float_t,dim,dim>;
1236  Eigen::LLT<Mat> lltX(X);
1237  Eigen::LLT<Mat> lltPsi(Psi);
1238  if((nu < dim) | (lltX.info() == Eigen::NumericalIssue) | (lltPsi.info() == Eigen::NumericalIssue)) return log ? -std::numeric_limits<float_t>::infinity() : 0.0;
1239  // https://stackoverflow.com/questions/35227131/eigen-check-if-matrix-is-positive-semi-definite
1240 
1241  float_t ldx (0.0); // log determinant of X
1242  float_t ldPsi (0.0);
1243  Mat Lx = lltX.matrixL(); // the lower diagonal L such that X = LL^T
1244  Mat Lpsi = lltPsi.matrixL();
1245  float_t logGammaNuOver2 = .25*dim*(dim-1)*log_pi<float_t>; // existence guaranteed when n > dim-1
1246 
1247  // add up log of diagonals of each Cholesky L
1248  for(size_t i = 0; i < dim; ++i){
1249  ldx += std::log(Lx(i,i));
1250  ldPsi += std::log(Lpsi(i,i));
1251  logGammaNuOver2 += std::lgamma(.5*(nu-i)); // recall j = i+1
1252  }
1253  ldx *= 2.0; // X = LL^T
1254  ldPsi *= 2.0;
1255 
1256  // TODO: this will probably be faster if you find an analogue...
1257  //float_t logDens = .5*nu*ldPsi - .5*(nu + dim + 1.0)*ldx - .5*(X.solve(Psi)).trace() - .5*nu*dim*std::log(2.0) - logGammaNuOver2;
1258  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;
1259 
1260 
1261  if(log){
1262  return logDens;
1263  }else{
1264  return std::exp(logDens);
1265  }
1266 }
1267 
1268 
1269 } //namespace rveval
1270 
1271 } //namespace pf
1272 
1273 #endif //RV_EVAL_H