coupling

coupling.c

#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);
}

coupling2.c

#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);
}