Skip to content

Instantly share code, notes, and snippets.

@dpiponi
Created May 1, 2021 19:46
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save dpiponi/dd73cd39e4e370cee585566e9dc98a13 to your computer and use it in GitHub Desktop.
coupling
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
int bernoulli(float p, float r)
{
return r > 1 - p;
}
int poisson(float lambda, float r)
{
float p = exp(-lambda);
for (int i = 0; ; ++i)
{
if ((r -= p) <= 0) return i;
p *= lambda / (i + 1);
}
}
/* 0 1 */
/* | 1-p | p | */
/* | ** *********************| */
/* | exp(-p) | pexp(-p) | p^2exp(-p)/2 | ... | */
int main()
{
srand48(time(0));
printf("%f\n", drand48());
#define N 8
float p[N] = { 0.1, 0.15, 0.15, 0.2, 0.2, 0.2 };
int t1 = 0, t2 = 0;
for (int i = 0; i < N; ++i)
{
float r = drand48();
t1 += bernoulli(p[i], r);
t2 += poisson(p[i], r);
}
printf("%d %d\n", t1, t2);
}
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
int bernoulli(float p, float r)
{
return r > 1 - p;
}
int poisson(float lambda, float r)
{
float p = exp(-lambda);
for (int i = 0; ; ++i)
{
if ((r -= p) <= 0) return i;
p *= lambda / (i + 1);
}
}
/* 0 1 */
/* | 1-p | p | */
/* | ** *********************| */
/* | exp(-p) | pexp(-p) | p^2exp(-p)/2 | ... | */
int main()
{
srand48(time(0));
printf("%f\n", drand48());
#define N 8
float p[N] = { 0.1, 0.15, 0.15, 0.2, 0.2, 0.2 };
int t1 = 0, t2 = 0;
for (int i = 0; i < N; ++i)
{
float r = drand48();
t1 += bernoulli(p[i], r);
t2 += poisson(p[i], r);
}
printf("%d %d\n", t1, t2);
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment