fork download
copy 
  1. #include <iostream>
  2. #include <cmath>
  3. #include <limits>
  4.  
  5.  
  6.  
  7. const double log_pi (1.1447298858494);
  8. const double log_two_pi (1.83787706640935);
  9. const double inv_sqrt_2pi(0.3989422804014327);
  10.  
  11. template<typename T = double>
  12. auto normal_distribution_density(T x, T mu = 0, T sigma = 1)
  13. {
  14.     auto const y = (x - mu) / sigma;
  15.     return 1 / (sigma * std::sqrt(2 * std::exp(log_pi))) * std::exp(-y * y / 2);
  16. }
  17.  
  18. // https://g...content-available-to-author-only...b.com/tbrown122387/pf/blob/master/include/rv_eval.h#L140template<typename T = double>
  19. template<typename float_t>
  20. float_t evalUnivNorm(float_t x, float_t mu, float_t sigma, bool log)
  21. {
  22.     float_t exponent = -.5*(x - mu)*(x-mu)/(sigma*sigma);
  23.     if( sigma > 0.0){
  24.         if(log){
  25.             return -std::log(sigma) - .5*log_two_pi + exponent;
  26.         }else{
  27.             return inv_sqrt_2pi * std::exp(exponent) / sigma;
  28.         }
  29.     }else{
  30.         if(log){
  31.             return -std::numeric_limits<float_t>::infinity();
  32.         }else{
  33.             return 0.0;
  34.         }
  35.     }
  36. }
  37.  
  38. template<typename T = double>
  39. auto wrapped_normal_distribution_density(T x, T mu = 0, T sigma = 1, T epsilon = 0)
  40. {
  41. 	T s = normal_distribution_density(x, mu, sigma);
  42. 	for (int k = 1;; ++k)
  43. 	{
  44. 		T a = normal_distribution_density(x + k, mu, sigma),
  45. 			b = normal_distribution_density(x - k, mu, sigma);
  46. 		s += a + b;
  47. 		if (a <= epsilon && b <= epsilon)
  48. 			break;
  49. 	}
  50. 	return s;
  51. }
  52.  
  53.  
  54. // https://g...content-available-to-author-only...b.com/tbrown122387/pf/blob/23d0e94ff8bd987997efafad66c120c07aaeddee/include/rv_eval.h#L120
  55. template<typename float_t>
  56. float_t log_sum_exp(float_t a, float_t b)
  57. {
  58.     float_t m = std::max(a,b);
  59.     return m + std::log(std::exp(a-m) + std::exp(b-m));
  60. }
  61.  
  62. template<typename T = double>
  63. auto wrapped_normal_distribution_density_taylor(T x, T mu = 0, T sigma = 1, T epsilon = 0)
  64. {
  65. 	// k=0
  66. 	T log_result = evalUnivNorm(x, mu, sigma, true);
  67. 	T last_iter_log_r;
  68. 	for (int absk = 1; absk < 1000; absk++)
  69. 	{
  70. 		last_iter_log_r = log_result;
  71. 		std::cout << "ay\n";
  72. 		log_result = log_sum_exp<T>(log_result, evalUnivNorm<T>(x + absk, mu, sigma, true));
  73. 		log_result = log_sum_exp<T>(log_result, evalUnivNorm<T>(x - absk, mu, sigma, true));
  74. 		if (last_iter_log_r == log_result)
  75. 			return std::exp(log_result);
  76. 	}
  77. }
  78.  
  79.  
  80.  
  81.  
  82. int main() {
  83.  
  84. 	std::cout << "your eval: " << wrapped_normal_distribution_density<double>(10.0) << "\n";
  85. 	std::cout << "my eval: " << wrapped_normal_distribution_density_taylor<double>(10.0) << "\n";
  86. 	return 0;
  87. }
Success #stdin #stdout 0s 4440KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
your eval: 1
my eval: ay
ay
ay
ay
ay
ay
ay
ay
ay
ay
ay
ay
ay
ay
ay
ay
ay
ay
ay
1
https://ideone.com/B2VLv3
language:
C++ (gcc 8.3)
created:
5 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!