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

Close

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

Download this file

KS.cc    79 lines (70 with data), 1.5 kB

 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 <config.h>
#include "KS.h"
#include <rng/RNG.h>
#include <cmath>
#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;
}
}
}