blob: 95da15ec9114b1e5849c1c72fefb1aa75bb18e66 (
plain)
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
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
|
#ifndef POISSONIER_H
#define POISSONIER_H
#include <random> /* for seed_seq */
#include <iostream>
#include <time.h> /* for clock_gettime */
#include <vector>
#ifdef BOOST_CPD
# include <boost/math/distributions/poisson.hpp>
#endif
static char const printable[64+1] = "ABCDEFGHIJKLMNOPQRSTUVWXYZ"
"abcdefghijklmnopqrstuvwxyz"
"0123456789._";
std::string base64(unsigned long int arg) {
std::string str;
for (;;) {
str = printable[arg & 63] + str;
arg >>= 6;
if (! arg) break;
}
return str;
}
std::string randomstr(unsigned int const nn=15) {
using namespace std;
// FIXME: ??? Could maybe use std::bitset here:
random_device rd;
string str;
while (str.length() < nn) {
unsigned int foo(rd());
int nch = (sizeof(foo) * 8) / 6;
for (int ii = 0; ii < nch; ii++) {
str += printable[foo & 63];
foo >>= 6;
}
}
return str;
}
class poissonier {
public:
std::mt19937 basegen;
std::string graine; // last seed fed to basegen
// constructor
poissonier () {
timespec time;
clock_gettime(CLOCK_REALTIME, &time);
graine = base64(time.tv_sec);
graine += base64(time.tv_nsec);
seed(graine);
}
double sample(double const lambda) {
// There is no way to change the intensity after the
// distro is constructed, so construct a new one each time:
return std::poisson_distribution<>(lambda)(basegen);
}
void seed(std::string const seed) {
// Need a temporary variable here, because the seed() function
// modifies it:
graine = seed;
std::seed_seq seed1(seed.begin(), seed.end());
basegen.seed(seed1);
}
std::string randomize() {
graine = randomstr(15);
seed(graine);
return graine;
}
void check(double const rate) {
using namespace std;
cout << "Checking: " << rate
<< " gives: " << sample(rate)
<< endl;
}
// Note this is stateless. Static member function.
//
// Log of the probability mass distribution:
static double lpmd(double const lambda, int const k){
// https://en.wikipedia.org/wiki/Poisson_distribution#Definition
if (lambda < 0) throw std::invalid_argument("negative Poisson intensity");
if (lambda == 0) {
if (k == 0) return 0; // 100% probability
return -INFINITY; // 0% probability
}
return k * log(lambda) - lambda - lgamma(k + 1.0);
}
// pmd = probability mass distribution
inline static double pmd(double const lambda, int const k){
return exp(lpmd(lambda, k));
}
#ifdef BOOST_CPD
// not implemented
// not needed, and would require boost
// cpd = cumulative probability distribution
static double cpd(double const lambda, int const k){
boost::math::poisson_distribution<> pd(lambda);
return boost::math::cpd(pd, k);
}
#endif
};
#endif /* POISSONIER_H */
|
