Learn how easy it is to sync an existing GitHub or Google Code repo to a SourceForge project!

## [01a998]: src / modules / glm / samplers / KS.cc Maximize Restore History

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78``` ```#include #include "KS.h" #include #include #define PI 3.141592653589793238462643383280 #define PI_SQUARE 9.86960440108936 using std::log; using std::exp; using std::pow; using std::sqrt; using std::fabs; static bool r_intvl(double u, double lambda) { double z = 1; double x = exp(-lambda/2); int j = 0; while(true) { j++; int j2 = (j+1)*(j+1); z -= j2 * pow(x, j2 - 1); if (z > u) return true; j++; j2 = (j+1)*(j+1); z += j2 * pow(x, j2 - 1); if (z <= u) return false; } } static bool l_intvl(double u, double lambda) { if (u==0) { return false; } double h = 0.5*log(2.0) + 2.5*log(PI) - 2.5*log(lambda) - PI_SQUARE/(2 * lambda) + 0.5*lambda; double logu = log(u); double z = 1; double x = exp(-PI_SQUARE/(2*lambda)); double k = lambda / PI_SQUARE; int j = 0; while(true) { j++; z -= k*pow(x, j*j-1); if (h + log(z) > logu) return true; j++; int j2 = (j+1)*(j+1); z += j2*pow(x,j2-1); if (h + log(z) <= logu) return false; } } namespace glm { double sample_lambda(double delta, RNG *rng) { delta = fabs(delta); while (true) { double y = rng->normal(); y = y * y; y = 1 + (y - sqrt(y*(4*delta + y)))/(2*delta); double u = rng->uniform(); double lambda = (u <= 1/(1+y)) ? delta/y : delta*y; u = rng->uniform(); bool ok = (lambda > 2.25) ? r_intvl(u, lambda) : l_intvl(u, lambda); if (ok) return lambda; } } } ```