pf
rv_eval.h
1 #ifndef RV_EVAL_H
2 #define RV_EVAL_H
3 
4 #include <cstddef> // std::size_t
5 #include <Eigen/Dense>
6 #include <iostream> // cerr
7 #include "boost/math/special_functions.hpp"
8 
9 
10 namespace rveval{
11 
15 
17 template<class T>
18 constexpr T inv_sqrt_2pi = T(0.3989422804014327);
19 
21 template<class T>
22 constexpr T sqrt_two_over_pi(0.797884560802865);
23 
25 template<class T>
26 constexpr T log_two_pi (1.83787706640935);
27 
29 template<class T>
30 constexpr T log_two_over_pi (-0.451582705289455);
31 
33 template<class T>
34 constexpr T log_pi (1.1447298858494);
35 
36 
40 
46 template<typename float_t>
47 float_t twiceFisher(float_t phi)
48 {
49  if ( (phi <= -1.0) || (phi >= 1.0) )
50  throw std::invalid_argument( "error: phi was not between -1 and 1" );
51  else
52  return std::log(1.0 + phi) - std::log(1.0 - phi);
53 }
54 
55 
56 
62 template<typename float_t>
63 float_t invTwiceFisher(float_t psi)
64 {
65  float_t ans = (1.0 - std::exp(psi)) / ( -1.0 - std::exp(psi) );
66 
67  if ( (ans <= -1.0) || (ans >= 1.0) )
68  std::cerr << "error: there was probably overflow for exp(psi) \n";
69 
70  return ans;
71 }
72 
73 
79 template<typename float_t>
80 float_t logit(float_t p)
81 {
82  if ( (p <= 0.0) || (p >= 1.0))
83  std::cerr << "error: p was not between 0 and 1 \n";
84 
85  return std::log(p) - std::log(1.0 - p);
86 }
87 
88 
94 template<typename float_t>
95 float_t inv_logit(float_t r)
96 {
97  float_t ans = 1.0/( 1.0 + std::exp(-r) );
98 
99  if ( (ans <= 0.0) || (ans >= 1.0))
100  std::cerr << "error: there was probably underflow for exp(-r) \n";
101 
102  return ans;
103 }
104 
105 
111 template<typename float_t>
112 float_t log_inv_logit(float_t r)
113 {
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));
116 }
117 
118 
125 template<typename float_t>
126 float_t log_sum_exp(float_t a, float_t b)
127 {
128  float_t m = std::max(a,b);
129  return m + std::log(std::exp(a-m) + std::exp(b-m));
130 }
131 
132 
136 
145 template<typename float_t>
146 float_t evalUnivNorm(float_t x, float_t mu, float_t sigma, bool log)
147 {
148  float_t exponent = -.5*(x - mu)*(x-mu)/(sigma*sigma);
149  if( sigma > 0.0){
150  if(log){
151  return -std::log(sigma) - .5*log_two_pi<float_t> + exponent;
152  }else{
153  return inv_sqrt_2pi<float_t> * std::exp(exponent) / sigma;
154  }
155  }else{
156  if(log){
157  return -std::numeric_limits<float_t>::infinity();
158  }else{
159  return 0.0;
160  }
161  }
162 }
163 
164 
173 template<typename float_t>
174 float_t evalUnivNorm_unnorm(float_t x, float_t mu, float_t sigma, bool log)
175 {
176  float_t exponent = -.5*(x - mu)*(x-mu)/(sigma*sigma);
177  if( sigma > 0.0){
178  if(log){
179  return exponent;
180  }else{
181  return std::exp(exponent);
182  }
183  }else{
184  if(log){
185  return -std::numeric_limits<float_t>::infinity();
186  }else{
187  return 0.0;
188  }
189  }
190 }
191 
192 
198 template<typename float_t>
199 float_t evalUnivStdNormCDF(float_t x) // john cook code
200 {
201  // constants
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;
208 
209  // Save the sign of x
210  int sign = 1;
211  if (x < 0)
212  sign = -1;
213  float_t xt = std::fabs(x)/std::sqrt(2.0);
214 
215  // A&S formula 7.1.26
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);
218 
219  return 0.5*(1.0 + sign*y);
220 }
221 
222 
231 template<typename float_t>
232 float_t evalUnivBeta(float_t x, float_t alpha, float_t beta, bool log)
233 {
234  if( (x > 0.0) && (x < 1.0) && (alpha > 0.0) && (beta > 0.0) ){ // x in support and parameters acceptable
235  if(log){
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);
237  }else{
238  return pow(x, alpha-1.0) * pow(1.0-x, beta-1.0) * std::tgamma(alpha + beta) / ( std::tgamma(alpha) * std::tgamma(beta) );
239  }
240 
241  }else{ //not ( x in support and parameters acceptable )
242  if(log){
243  return -std::numeric_limits<float_t>::infinity();
244  }else{
245  return 0.0;
246  }
247  }
248 }
249 
250 
259 template<typename float_t>
260 float_t evalUnivBeta_unnorm(float_t x, float_t alpha, float_t beta, bool log)
261 {
262  if( (x > 0.0) && (x < 1.0) && (alpha > 0.0) && (beta > 0.0) ){ // x in support and parameters acceptable
263  if(log){
264  return (alpha - 1.0)*std::log(x) + (beta - 1.0) * std::log(1.0 - x);
265  }else{
266  return pow(x, alpha-1.0) * pow(1.0-x, beta-1.0);
267  }
268 
269  }else{ //not ( x in support and parameters acceptable )
270  if(log){
271  return -std::numeric_limits<float_t>::infinity();
272  }else{
273  return 0.0;
274  }
275  }
276 }
277 
278 
287 template<typename float_t>
288 float_t evalUnivInvGamma(float_t x, float_t alpha, float_t beta, bool log)
289 {
290  if ( (x > 0.0) && (alpha > 0.0) && (beta > 0.0) ){ // x in support and acceptable parameters
291  if (log){
292  return alpha * std::log(beta) - std::lgamma(alpha) - (alpha + 1.0)*std::log(x) - beta/x;
293  }else{
294  return pow(x, -alpha-1.0) * exp(-beta/x) * pow(beta, alpha) / std::tgamma(alpha);
295  }
296  }else{ // not ( x in support and acceptable parameters )
297  if (log){
298  return -std::numeric_limits<float_t>::infinity();
299  }else{
300  return 0.0;
301  }
302  }
303 }
304 
305 
314 template<typename float_t>
315 float_t evalUnivInvGamma_unnorm(float_t x, float_t alpha, float_t beta, bool log)
316 {
317  if ( (x > 0.0) && (alpha > 0.0) && (beta > 0.0) ){ // x in support and acceptable parameters
318  if (log){
319  return (-alpha - 1.0)*std::log(x) - beta/x;
320  }else{
321  return pow(x, -alpha-1.0) * exp(-beta/x);
322  }
323  }else{ // not ( x in support and acceptable parameters )
324  if (log){
325  return -std::numeric_limits<float_t>::infinity();
326  }else{
327  return 0.0;
328  }
329  }
330 }
331 
332 
340 template<typename float_t>
341 float_t evalUnivHalfNorm(float_t x, float_t sigmaSqd, bool log)
342 {
343  if( (x >= 0.0) && (sigmaSqd > 0.0)){
344  if (log){
345  return .5*log_two_over_pi<float_t> - .5*std::log(sigmaSqd) - .5*x*x / sigmaSqd;
346  }else{
347  return std::exp(-.5*x*x/sigmaSqd) * sqrt_two_over_pi<float_t> / std::sqrt(sigmaSqd);
348  }
349  }else{
350  if (log){
351  return -std::numeric_limits<float_t>::infinity();
352  }else{
353  return 0.0;
354  }
355  }
356 }
357 
358 
366 template<typename float_t>
367 float_t evalUnivHalfNorm_unnorm(float_t x, float_t sigmaSqd, bool log)
368 {
369  if( (x >= 0.0) && (sigmaSqd > 0.0)){
370  if (log){
371  return -.5*x*x / sigmaSqd;
372  }else{
373  return std::exp(-.5*x*x/sigmaSqd);
374  }
375  }else{
376  if (log){
377  return -std::numeric_limits<float_t>::infinity();
378  }else{
379  return 0.0;
380  }
381  }
382 }
383 
384 
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)
397 {
398  if( (sigma > 0.0) && (lower <= x) & (x <= upper) ){
399  if(log){
400  return evalUnivNorm(x, mu, sigma, true)
401  - std::log( evalUnivStdNormCDF((upper-mu)/sigma) - evalUnivStdNormCDF((lower-mu)/sigma));
402  }else{
403  return evalUnivNorm(x,mu,sigma,false)
404  / ( evalUnivStdNormCDF((upper-mu)/sigma) - evalUnivStdNormCDF((lower-mu)/sigma) );
405  }
406 
407  }else{
408  if (log){
409  return -std::numeric_limits<float_t>::infinity();
410  }else{
411  return 0.0;
412  }
413  }
414 }
415 
416 
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)
429 {
430  if( (sigma > 0.0) && (lower <= x) & (x <= upper) ){
431  if(log){
432  return evalUnivNorm_unnorm(x, mu, sigma, true);
433  }else{
434  return evalUnivNorm_unnorm(x, mu, sigma, false);
435  }
436  }else{
437  if (log){
438  return -std::numeric_limits<float_t>::infinity();
439  }else{
440  return 0.0;
441  }
442  }
443 }
444 
445 
454 template<typename float_t>
455 float_t evalLogitNormal(float_t x, float_t mu, float_t sigma, bool log)
456 {
457  if( (x >= 0.0) && (x <= 1.0) && (sigma > 0.0)){
458 
459  float_t exponent = -.5*(logit(x) - mu)*(logit(x) - mu) / (sigma*sigma);
460  if(log){
461  return -std::log(sigma) - .5*log_two_pi<float_t> - std::log(x) - std::log(1.0-x) + exponent;
462  }else{
463  return inv_sqrt_2pi<float_t> * std::exp(exponent) / (x * (1.0-x) * sigma);
464  }
465  }else{
466  if(log){
467  return -std::numeric_limits<float_t>::infinity();
468  }else{
469  return 0.0;
470  }
471  }
472 }
473 
474 
483 template<typename float_t>
484 float_t evalLogitNormal_unnorm(float_t x, float_t mu, float_t sigma, bool log)
485 {
486  if( (x >= 0.0) && (x <= 1.0) && (sigma > 0.0)){
487 
488  float_t exponent = -.5*(logit(x) - mu)*(logit(x) - mu) / (sigma*sigma);
489  if(log){
490  return -std::log(x) - std::log(1.0-x) + exponent;
491  }else{
492  return std::exp(exponent) / x / (1.0-x);
493  }
494  }else{
495  if(log){
496  return -std::numeric_limits<float_t>::infinity();
497  }else{
498  return 0.0;
499  }
500  }
501 }
502 
503 
512 template<typename float_t>
513 float_t evalTwiceFisherNormal(float_t x, float_t mu, float_t sigma, bool log)
514 {
515 
516  // https://stats.stackexchange.com/questions/321905/what-is-the-name-of-this-random-variable/321907#321907
517  if( (x >= -1.0) && (x <= 1.0) && (sigma > 0.0)){
518 
519  float_t exponent = std::log((1.0+x)/(1.0-x)) - mu;
520  exponent = -.5*exponent*exponent/sigma/sigma;
521  if(log){
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;
523  }else{
524  return inv_sqrt_2pi<float_t> * 2.0 * std::exp(exponent)/( (1.0-x)*(1.0+x)*sigma );
525  }
526  }else{
527  if(log){
528  return -std::numeric_limits<float_t>::infinity();
529  }else{
530  return 0.0;
531  }
532  }
533 }
534 
535 
544 template<typename float_t>
545 float_t evalTwiceFisherNormal_unnorm(float_t x, float_t mu, float_t sigma, bool log)
546 {
547 
548  // https://stats.stackexchange.com/questions/321905/what-is-the-name-of-this-random-variable/321907#321907
549  if( (x >= -1.0) && (x <= 1.0) && (sigma > 0.0)){
550 
551  float_t exponent = std::log((1.0+x)/(1.0-x)) - mu;
552  exponent = -.5*exponent*exponent/sigma/sigma;
553  if(log){
554  return -std::log(1.0+x) - std::log(1.0-x) + exponent;
555  }else{
556  return std::exp(exponent)/(1.0-x)/(1.0+x);
557  }
558  }else{
559  if(log){
560  return -std::numeric_limits<float_t>::infinity();
561  }else{
562  return 0.0;
563  }
564  }
565 }
566 
567 
576 template<typename float_t>
577 float_t evalLogNormal(float_t x, float_t mu, float_t sigma, bool log)
578 {
579  if( (x > 0.0) && (sigma > 0.0)){
580 
581  float_t exponent = std::log(x)-mu;
582  exponent = -.5 * exponent * exponent / sigma / sigma;
583  if(log){
584  return -std::log(x) - std::log(sigma) - .5*log_two_pi<float_t> + exponent;
585  }else{
586  return inv_sqrt_2pi<float_t> * std::exp(exponent)/(sigma*x);
587  }
588  }else{
589  if(log){
590  return -std::numeric_limits<float_t>::infinity();
591  }else{
592  return 0.0;
593  }
594  }
595 }
596 
597 
606 template<typename float_t>
607 float_t evalLogNormal_unnorm(float_t x, float_t mu, float_t sigma, bool log)
608 {
609  if( (x > 0.0) && (sigma > 0.0)){
610 
611  float_t exponent = std::log(x)-mu;
612  exponent = -.5 * exponent * exponent / sigma / sigma;
613  if(log){
614  return -std::log(x) + exponent;
615  }else{
616  return std::exp(exponent)/x;
617  }
618  }else{
619  if(log){
620  return -std::numeric_limits<float_t>::infinity();
621  }else{
622  return 0.0;
623  }
624  }
625 }
626 
627 
636 template<typename float_t>
637 float_t evalUniform(float_t x, float_t lower, float_t upper, bool log)
638 {
639 
640  if( (x > lower) && (x <= upper)){
641 
642  float_t width = upper-lower;
643  if(log){
644  return -std::log(width);
645  }else{
646  return 1.0/width;
647  }
648  }else{
649  if(log){
650  return -std::numeric_limits<float_t>::infinity();
651  }else{
652  return 0.0;
653  }
654  }
655 
656 }
657 
658 
667 template<typename float_t>
668 float_t evalUniform_unnorm(float_t x, float_t lower, float_t upper, bool log)
669 {
670 
671  if( (x > lower) && (x <= upper)){
672 
673  if(log){
674  return 0.0;
675  }else{
676  return 1.0;
677  }
678  }else{
679  if(log){
680  return -std::numeric_limits<float_t>::infinity();
681  }else{
682  return 0.0;
683  }
684  }
685 
686 }
687 
688 
698 template<typename float_t>
699 float_t evalScaledT(float_t x, float_t mu, float_t sigma, float_t dof, bool log)
700 {
701 
702  if( (sigma > 0.0) && (dof > 0.0) ){
703 
704  float_t zscore = (x-mu)/sigma;
705  float_t lmt = - .5*(dof+1.0)*std::log(1.0 + (zscore*zscore)/dof);
706  if(log)
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;
708  else
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);
710  } else{
711  if(log)
712  return -std::numeric_limits<float_t>::infinity();
713  else
714  return 0.0;
715  }
716 }
717 
718 
728 template<typename float_t>
729 float_t evalScaledT_unnorm(float_t x, float_t mu, float_t sigma, float_t dof, bool log)
730 {
731  if( (sigma > 0.0) && (dof > 0.0) ){
732 
733  float_t zscore = (x-mu)/sigma;
734  float_t lmt = - .5*(dof+1.0)*std::log(1.0 + (zscore*zscore)/dof);
735  if(log)
736  return lmt;
737  else
738  return std::exp(lmt);
739  } else{
740  if(log)
741  return -std::numeric_limits<float_t>::infinity();
742  else
743  return 0.0;
744  }
745 }
746 
747 
755 template<typename int_t, typename float_t>
756 float_t evalDiscreteUnif(int_t x, int k, bool log)
757 {
758  if( (1 <= x) && (x <= k) ){
759  if(log){
760  return -std::log(static_cast<float_t>(k));
761  }else{
762  return 1.0 / static_cast<float_t>(k);
763  }
764  }else{ // x not in support
765  if(log){
766  return -std::numeric_limits<float_t>::infinity();
767  }else{
768  return 0.0;
769  }
770  }
771 }
772 
773 
781 template<typename int_t, typename float_t>
782 float_t evalDiscreteUnif_unnorm(int_t x, int_t k, bool log)
783 {
784  if( (1 <= x) && (x <= k) ){
785  if(log){
786  return 0.0;
787  }else{
788  return 1.0;
789  }
790  }else{
791  if(log){
792  return -std::numeric_limits<float_t>::infinity();
793  }else{
794  return 0.0;
795  }
796  }
797 }
798 
799 
806 template<typename int_t, typename float_t>
807 float_t evalBernoulli(int_t x, float_t p, bool log)
808 {
809  if( ((x == 0) || (x == 1)) && ( (0.0 <= p) && (p <= 1.0) ) ){ // if valid x and valid p
810  if(log){
811  return (x==1) ? std::log(p) : std::log(1.0-p);
812  }else{
813  return (x==1) ? p : (1.0-p);
814  }
815  }else{ // either invalid x or invalid p
816  if(log){
817  return -std::numeric_limits<float_t>::infinity();
818  }else{
819  return 0.0;
820  }
821  }
822 }
823 
824 
833 template<typename int_t, typename float_t>
834 float_t evalSkellam(int_t x, float_t mu1, float_t mu2, bool log)
835 {
836  if( (mu1 > 0) && (mu2 > 0) ){
837 
838  // much of this function is adapted from
839  // https://github.com/stan-dev/math/blob/9b2e93ba58fa00521275b22a190468ab22f744a3/stan/math/prim/fun/log_modified_bessel_first_kind.hpp
840 
841  // step 1: calculate log I_k(2\sqrt{mu_1 mu_2}) using log_sum_exp
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());
845 
846  if (x == 0) {
847  // modified from Boost's bessel_i0_imp in the double precision case,
848  // which refers to:
849  // Modified Bessel function of the first kind of order zero
850  // we use the approximating forms derived in:
851  // "Rational Approximations for the Modified Bessel Function of the
852  // First Kind -- I0(x) for Computations with Double Precision"
853  // by Pavel Holoborodko, see
854  // http://www.advanpix.com/2015/11/11/rational-approximations-for-the-modified-bessel-function-of-the-first-kind-i0-computations-double-precision
855  // The actual coefficients used are [Boost's] own, and extend
856  // Pavel's work to precisions other than double.
857 
858  if (z < 7.75) {
859  // Bessel I0 over[10 ^ -16, 7.75]
860  // Max error in interpolated form : 3.042e-18
861  // Max Error found at double precision = Poly : 5.106609e-16
862  // Cheb : 5.239199e-16
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)));
873 
874  }else if (z < 500) {
875  // Max error in interpolated form : 1.685e-16
876  // Max Error found at double precision = Poly : 2.575063e-16
877  // Cheb : 2.247615e+00
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};
890 
891  log_I = z + std::log(evaluate_polynomial(P, 1.0/z)) - 0.5*std::log(z);
892 
893  }else{
894 
895  // Max error in interpolated form : 2.437e-18
896  // Max Error found at double precision = Poly : 1.216719e-16
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);
901 
902  }
903 
904  }else if(std::abs(x)==1){
905 
906  // modified from Boost's bessel_i1_imp in the double precision case
907  // see credits above in the v == 0 case
908  if (z < 7.75) {
909  // Bessel I0 over[10 ^ -16, 7.75]
910  // Max error in interpolated form: 5.639e-17
911  // Max Error found at double precision = Poly: 1.795559e-16
912 
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};
921 
922  float_t a = mu1*mu2;
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);
925 
926  }else if (z < 500) {
927  // Max error in interpolated form: 1.796e-16
928  // Max Error found at double precision = Poly: 2.898731e-16
929 
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);
943  }else {
944 
945  // Max error in interpolated form: 1.320e-19
946  // Max Error found at double precision = Poly: 7.065357e-17
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);
952  }
953 
954  }else if( z > 100){
955 
956  // Boost does something like this in asymptotic_bessel_i_large_x
957  float_t lim = std::pow((x*x + 2.5) / (2 * z), 3) / 24;
958 
959  if (lim < std::numeric_limits<float_t>::epsilon()* 10) {
960  float_t s = 1;
961  float_t mu = 4*x*x;
962  float_t ex = 8 * z;
963  float_t num = mu - 1;
964  float_t denom = ex;
965  s -= num / denom;
966  num *= mu - 9;
967  denom *= ex * 2;
968  s += num / denom;
969  num *= mu - 25;
970  denom *= ex * 3;
971  s -= num / denom;
972  log_I = z - .5*std::log(z) - .5*log_two_pi<float_t> + std::log(s);
973  }
974  }else{
975 
976  // just do the sum straightforwardly
977  // m=0
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;
984 
985  // m > 0
986  float_t m = 1.0;
987  float_t last_iter_log_I;
988  do{
989  first += lm1m2;
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);
994  m++;
995  if(m > 1000)
996  std::cout << "first, second, third, mu1, mu2: " << first << ", " << second << ", " << third << ", " << mu1 << ", " << mu2 << "\n";
997  }while(log_I != last_iter_log_I);
998 
999  }
1000 
1001  // step 2: add the easy parts to get the overall pmf evaluation
1002  float_t log_mass = -mu1 - mu2 + .5*x*(std::log(mu1) - std::log(mu2)) + log_I;
1003  //std::cout << "guaranteed: " << std::log(boost::math::cyl_bessel_i<int_t, float_t>(x,2.0*std::sqrt(mu1*mu2)))
1004  // << "\nmine: " << log_I << "\n";
1005 
1006  // step 3: handle log/nonlog particulars
1007  if(log) {
1008  return log_mass;
1009  }else {
1010  return std::exp(log_mass);
1011  }
1012  }else {
1013  if(log) {
1014  return -std::numeric_limits<float_t>::infinity();
1015  }else {
1016  return 0.0;
1017  }
1018  }
1019 }
1020 
1024 
1025 
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,
1042  bool log = false)
1043 {
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; // if not pd return 0 dens
1047  Mat L = lltM.matrixL(); // the lower diagonal L such that M = LL^T
1048  float_t quadform = (lltM.solve(x-meanVec)).squaredNorm();
1049  float_t ld (0.0); // calculate log-determinant using cholesky decomposition too
1050  // add up log of diagnols of Cholesky L
1051  for(size_t i = 0; i < dim; ++i){
1052  ld += std::log(L(i,i));
1053  }
1054  ld *= 2; // covMat = LL^T
1055 
1056  float_t logDens = -.5*log_two_pi<float_t> * dim - .5*ld - .5*quadform;
1057 
1058  if(log){
1059  return logDens;
1060  }else{
1061  return std::exp(logDens);
1062  }
1063 }
1064 
1065 
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,
1082  const float_t& dof,
1083  bool log = false)
1084 {
1085  if(dof <= 0.0) return log ? -std::numeric_limits<float_t>::infinity() : 0.0; // degrees of freedom must be positive
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; // if not pd return 0 dens
1089  Mat L = lltM.matrixL(); // the lower diagonal L such that M = LL^T
1090  float_t quadform = (lltM.solve(x-locVec)).squaredNorm();
1091  float_t ld (0.0); // calculate log-determinant using cholesky decomposition too
1092  // add up log of diagnols of Cholesky L
1093  for(size_t i = 0; i < dim; ++i){
1094  ld += std::log(L(i,i));
1095  }
1096  ld *= 2; // shapeMat = LL^T
1097 
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 );
1100 
1101  if(log){
1102  return logDens;
1103  }else{
1104  return std::exp(logDens);
1105  }
1106 }
1107 
1108 
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,
1126  bool log = false)
1127 {
1128 
1129  using bigmat = Eigen::Matrix<float_t,bigd,bigd>;
1130  using smallmat = Eigen::Matrix<float_t,smalld,smalld>;
1131 
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(); // LL' = Sig^{-1}
1138  float_t quadform = (L * (x-meanVec)).squaredNorm();
1139  if (log){
1140 
1141  // calculate log-determinant using cholesky decomposition (assumes symmetric and positive definite)
1142  float_t halfld (0.0);
1143  // add up log of diagnols of Cholesky L
1144  for(size_t i = 0; i < bigd; ++i){
1145  halfld += std::log(L(i,i));
1146  }
1147 
1148  return -.5*log_two_pi<float_t> * bigd + halfld - .5*quadform;
1149 
1150 
1151  }else{ // not the log density
1152  float_t normConst = std::pow(inv_sqrt_2pi<float_t>, bigd) * L.determinant();
1153  return normConst * std::exp(-.5* quadform);
1154  }
1155 
1156 }
1157 
1158 
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,
1175  bool log = false)
1176 {
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;
1181  // https://stackoverflow.com/questions/35227131/eigen-check-if-matrix-is-positive-semi-definite
1182 
1183  float_t ldx (0.0); // log determinant of X
1184  float_t ldvinv (0.0);
1185  Mat Lx = lltX.matrixL(); // the lower diagonal L such that X = LL^T
1186  Mat Lvi = lltVinv.matrixL();
1187  float_t logGammaNOver2 = .25*dim*(dim-1)*log_pi<float_t>; // existence guaranteed when n > dim-1
1188 
1189  // add up log of diagonals of each Cholesky L
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)); // recall j = i+1
1194  }
1195  ldx *= 2.0; // X = LL^T
1196  ldvinv *= 2.0;
1197 
1198  float_t logDens = .5*(n - dim -1)*ldx - .5*(Vinv*X).trace() - .5*n*dim*std::log(2.0) + .5*n*ldvinv - logGammaNOver2;
1199 
1200  if(log){
1201  return logDens;
1202  }else{
1203  return std::exp(logDens);
1204  }
1205 }
1206 
1207 
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,
1224  bool log = false)
1225 {
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;
1230  // https://stackoverflow.com/questions/35227131/eigen-check-if-matrix-is-positive-semi-definite
1231 
1232  float_t ldx (0.0); // log determinant of X
1233  float_t ldPsi (0.0);
1234  Mat Lx = lltX.matrixL(); // the lower diagonal L such that X = LL^T
1235  Mat Lpsi = lltPsi.matrixL();
1236  float_t logGammaNuOver2 = .25*dim*(dim-1)*log_pi<float_t>; // existence guaranteed when n > dim-1
1237 
1238  // add up log of diagonals of each Cholesky L
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)); // recall j = i+1
1243  }
1244  ldx *= 2.0; // X = LL^T
1245  ldPsi *= 2.0;
1246 
1247  // TODO: this will probably be faster if you find an analogue...
1248  //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;
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;
1250 
1251 
1252  if(log){
1253  return logDens;
1254  }else{
1255  return std::exp(logDens);
1256  }
1257 }
1258 
1259 
1260 } //namespace rveval
1261 
1262 
1263 #endif //RV_EVAL_H
Definition: rv_eval.h:10