Fleet  0.0.9
Inference in the LOT
Random.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include <fstream>
4 #include <random>
5 #include <functional>
6 
7 #include "Errors.h"
8 #include "Numerics.h"
9 #include "Rng.h"
10 
11 double uniform() {
16  std::uniform_real_distribution<double> d(0.0, 1.0);
17  return d(rng);
18 }
19 
20 double uniform(const double a, const double b) {
21  std::uniform_real_distribution<double> d(a,b);
22  return d(rng);
23 }
24 
25 bool flip(float p = 0.5) {
31  return uniform() < p;
32 }
33 
34 
35 double random_normal(double mu=0, double sd=1.0) {
36  std::normal_distribution<float> normal(0.0, 1.0);
37  return normal(rng)*sd + mu;
38 }
39 
40 template<typename T>
41 T normal_lpdf(T x, T mu=0.0, T sd=1.0) {
50  //https://stackoverflow.com/questions/10847007/using-the-gaussian-probability-density-function-in-c
51  const T linv_sqrt_2pi = -0.5*log(2*pi*sd*sd);
52  const T z = (x-mu)/sd;
53  return linv_sqrt_2pi - z*z / 2.0;
54 }
55 
56 template<typename T>
57 T normal_cdf(T x, T mu, T sd) {
58  T z = (x-mu)/sd;
59  return std::erfc(-z/std::sqrt(2))/2;
60 }
61 
62 template<typename T>
63 T random_t(T nu) {
64  std::student_t_distribution<double> dist(nu);
65  return dist(rng);
66 }
67 
68 template<typename T>
69 T t_lpdf(T x, T nu) {
70  return mylgamma((nu+1)/2) - log(nu*pi)/2 - mylgamma(nu/2) +
71  log(1+x*x/nu)*(-(nu+1)/2);
72 }
73 
74 
75 double random_cauchy() {
81  return tan(pi*(uniform()-0.5));
82 }
83 
84 
85 template<typename T>
86 T cauchy_lpdf(T x, T loc=0.0, T gamma=1.0) {
95  return -log(pi) - log(gamma) - log(1+((x-loc)/gamma)*((x-loc)/gamma));
96 }
97 
98 
99 template<typename T>
100 T laplace_lpdf(T x, T mu=0.0, T b=1.0) {
109  //https://stackoverflow.com/questions/10847007/using-the-gaussian-probability-density-function-in-c
110  return -log(2*b) - std::abs(x-mu)/b;
111 }
112 
113 
114 template<typename T>
115 T random_exponential(T l=1.0) {
123  return -log(uniform())/l;
124 }
125 
126 
127 template<typename T>
128 T random_laplace(T mu=0.0, T b=1.0) {
136  return mu + random_exponential(b) - random_exponential(b);
137 }
138 
139 
140 double geometric_lpdf(size_t k, double p) {
141  return log(p) + (k-1)*log(1-p);
142 }
143 
144 
145 template<typename t>
146 std::vector<t> random_multinomial(t a, size_t len) {
154  std::gamma_distribution<t> g(a, 1.0);
155  std::vector<t> out(len);
156  t total = 0.0;
157  for(size_t i =0;i<len;i++){
158  out[i] = g(rng);
159  total += out[i];
160  }
161  for(size_t i=0;i<len;i++){
162  out[i] /= total;
163  }
164  return out;
165 }
166 
167 template<typename T>
168 T myrandom(T max) {
174  assert(max>0);
175  std::uniform_int_distribution<T> r(0, max-1);
176  return r(rng);
177 }
178 
179 template<typename T>
180 T myrandom(T min, T max) {
188  std::uniform_int_distribution<T> r(min,max-1);
189  return r(rng);
190 }
191 
192 
193 template<typename t, typename T>
194 double sample_z(const T& s, const std::function<double(const t&)>& f) {
202  double z = 0.0;
203  for(auto& x: s) {
204  auto fx = f(x);
205  if(not std::isnan(fx))
206  z += f(x);
207  }
208  return z;
209 }
210 
211 template<typename t, typename T>
212 std::pair<t*,double> sample(const T& s, double z, const std::function<double(const t&)>& f = [](const t& v){return 1.0;}) {
213  // this takes a collection T of elements t, and a function f
214  // which assigns them each a probability, and samples from them according
215  // to the probabilities in f. The probability that is returned is only the probability
216  // of selecting that *element* (index), not the probability of selecting anything equal to it
217  // (i.e. we defaultly don't double-count equal options). For that, use p_sample_eq below
218  // here z is the normalizer, z = sum_x f(x)
219  assert(z > 0 && "*** Cannot call sample with zero normalizer -- is s empty?");
220  double r = z * uniform();
221  for(auto& x : s) {
222 
223  auto fx = f(x);
224  if(std::isnan(fx)) continue; // treat as zero prob
225  assert(fx >= 0.0);
226  r -= fx;
227  if(r <= 0.0)
228  return std::make_pair(const_cast<t*>(&x), log(fx)-log(z));
229  }
230 
231  throw YouShouldNotBeHereError("*** Should not get here in sampling");
232 }
233 
234 template<typename t, typename T>
235 std::pair<t*,double> sample(const T& s, const std::function<double(const t&)>& f = [](const t& v){return 1.0;}) {
236  // An interface to sample that computes the normalizer for you
237  return sample<t,T>(s, sample_z<t,T>(s,f), f);
238 }
239 
240 
241 std::pair<int,double> sample_int(unsigned int max, const std::function<double(const int)>& f = [](const int v){return 1.0;}) {
242  // special form for where the ints (e.g. indices) Are what f takes)
243  double z = 0.0;
244  for(size_t i=0;i<max;i++){
245  z += f(i);
246  }
247 
248  double r = z * uniform();
249 
250  for(size_t i=0;i<max;i++){
251  double fx = f(i);
252  if(std::isnan(fx)) continue; // treat as zero prob
253  assert(fx >= 0.0);
254  r -= fx;
255  if(r <= 0.0)
256  return std::make_pair(i, log(fx)-log(z));
257  }
258 
259  throw YouShouldNotBeHereError("*** Should not get here in sampling");
260 }
261 
262 
263 std::pair<int,double> sample_int_lp(unsigned int max, const std::function<double(const int)>& f = [](const int v){return 1.0;}) {
264  // special form for where the ints (e.g. indices) Are what f takes)
265  double z = -infinity;
266  for(size_t i=0;i<max;i++){
267  double fx = f(i);
268  if(std::isnan(fx)) continue; // treat as zero prob
269  z = logplusexp(z, fx);
270  }
271 
272  double r = z + log(uniform());
273 
274  double fz = -infinity;
275  for(size_t i=0;i<max;i++){
276  double fx = f(i);
277  if(std::isnan(fx)) continue; // treat as zero prob
278  //assert(fx <= 0.0); // we actually can have lp > 0 if we wanted...
279  fz = logplusexp(fz, fx);
280  if(r <= fz)
281  return std::make_pair(i, fx-z);
282  }
283 
284  throw YouShouldNotBeHereError("*** Should not get here in sampling");
285 }
286 
287 
288 template<typename t, typename T>
289 std::pair<size_t,double> arg_max(const T& s, const std::function<double(const t&)>& f) {
290  // Same interface as sample but choosing the max
291  double mx = -infinity;
292  size_t max_idx = 0;
293  size_t idx = 0;
294  for(auto& x : s) {
295  double fx = f(x);
296  if(fx > mx) {
297  mx = fx;
298  max_idx = idx;
299  }
300  idx++;
301  }
302  return std::make_pair(max_idx, mx);
303 }
304 
305 
306 std::pair<size_t,double> arg_max_int(unsigned int max, const std::function<double(const int)>& f) {
307  double mx = -infinity;
308  size_t max_idx=0;
309  for(size_t i=0;i<max;i++){
310  double fx = f(i);
311  if(std::isnan(fx)) continue; // treat as zero prob
312  if(fx > mx) {
313  max_idx = i;
314  mx = fx;
315  }
316  }
317  return std::make_pair(max_idx, mx);
318 }
319 
320 template<typename t, typename T>
321 std::pair<t*,double> max_of(const T& s, const std::function<double(const t&)>& f) {
322  // Same interface as sample but choosing the max
323  double mx = -infinity;
324  t* out = nullptr;
325  for(auto& x : s) {
326  double fx = f(x);
327  if(fx > mx) {
328  mx = fx;
329  out = const_cast<t*>(&x);
330  }
331  }
332  return std::make_pair(out, mx);
333 }
334 
335 template<typename t, typename T>
336 std::pair<t*,double> sample(const T& s, double(*f)(const t&)) {
337  std::function sf = f;
338  return sample<t,T>(s,sf);
339 }
340 
341 template<typename t, typename T>
342 double lp_sample_eq(const t& x, const T& s, std::function<double(const t&)>& f = [](const t& v){return 1.0;}) {
343  // the probability of sampling *anything* equal to x out of s (including double counts)
344 
345  double z = sample_z(s,f);
346  double px = 0.0; // find the probability of anything *equal* to x
347  for(const auto& y : s) {
348  if(y == x)
349  px += f(y);
350  }
351 
352  assert(px <= z);
353 
354  return log(px)-log(z);
355 }
356 
357 template<typename t, typename T>
358 double lp_sample_eq(const t& x, const T& s, double(*f)(const t&)) {
359  std::function sf = f;
360  return lp_sample_eq<t,T>(x,s,sf);
361 }
362 
363 template<typename t, typename T>
364 double lp_sample_one(const t& x, const T& s, std::function<double(const t&)>& f = [](const t& v){return 1.0;}) {
365  // probability of sampling just this one
366  // NOTE: does not check that x is in s!
367 
368  double z = sample_z(s,f);
369  return log(f(x))-log(z);
370 }
371 template<typename t, typename T>
372 double lp_sample_one(const t& x, const T& s, double(*f)(const t&)) {
373  std::function sf = f;
374  return lp_sample_one<t,T>(x,s,sf);
375 }
T myrandom(T max)
Definition: Random.h:168
double lp_sample_one(const t &x, const T &s, std::function< double(const t &)> &f=[](const t &v){return 1.0;})
Definition: Random.h:364
T random_t(T nu)
Definition: Random.h:63
std::pair< size_t, double > arg_max(const T &s, const std::function< double(const t &)> &f)
Definition: Random.h:289
T t_lpdf(T x, T nu)
Definition: Random.h:69
std::pair< t *, double > sample(const T &s, double z, const std::function< double(const t &)> &f=[](const t &v){return 1.0;})
Definition: Random.h:212
Rng rng
double uniform()
Definition: Random.h:11
constexpr double pi
Definition: Numerics.h:20
std::pair< int, double > sample_int_lp(unsigned int max, const std::function< double(const int)> &f=[](const int v){return 1.0;})
Definition: Random.h:263
double mylgamma(double v)
Definition: Numerics.h:143
T laplace_lpdf(T x, T mu=0.0, T b=1.0)
Definition: Random.h:100
std::pair< int, double > sample_int(unsigned int max, const std::function< double(const int)> &f=[](const int v){return 1.0;})
Definition: Random.h:241
double geometric_lpdf(size_t k, double p)
Definition: Random.h:140
Definition: Errors.h:18
bool flip(float p=0.5)
Definition: Random.h:25
T random_exponential(T l=1.0)
Definition: Random.h:115
std::vector< t > random_multinomial(t a, size_t len)
Definition: Random.h:146
T random_laplace(T mu=0.0, T b=1.0)
Definition: Random.h:128
double random_normal(double mu=0, double sd=1.0)
Definition: Random.h:35
double random_cauchy()
Definition: Random.h:75
constexpr double infinity
Definition: Numerics.h:18
T logplusexp(const T a, const T b)
Definition: Numerics.h:89
double lp_sample_eq(const t &x, const T &s, std::function< double(const t &)> &f=[](const t &v){return 1.0;})
Definition: Random.h:342
std::pair< t *, double > max_of(const T &s, const std::function< double(const t &)> &f)
Definition: Random.h:321
T normal_cdf(T x, T mu, T sd)
Definition: Random.h:57
T cauchy_lpdf(T x, T loc=0.0, T gamma=1.0)
Definition: Random.h:86
std::pair< size_t, double > arg_max_int(unsigned int max, const std::function< double(const int)> &f)
Definition: Random.h:306
T normal_lpdf(T x, T mu=0.0, T sd=1.0)
Definition: Random.h:41
double sample_z(const T &s, const std::function< double(const t &)> &f)
Definition: Random.h:194