Appendix - Monte Carlo process for calculating confidence intervals.

If you check the textbooks, confidence intervals are a little ambiguously defined. Where they are defined in terms of a recipe, they usually refer to samples from a normal population or for a large number of samples. Clearly the situation here is neither.

Looking closely at the assumptions in the standard recipes shows two things :- They assume the sample size is large. They use the sample statistic as an estimator of the population statistic at some point.

One can get a handle on this situation by inverting the problem as follows :- Assume that the probability of finding the mouth open is p, about which we know nothing, (ie. p is a random variable which has a flat distribution between [0,1]), what is the probability distribution of p GIVEN that in n observations, r of them showed the mouth was open.

One can solve this by:-

Write a computer program to :-

1. Generate a random number p distributed uniformly between 0 and 1.
2. Generate n random (0,1) events with probability p of a 1
3. if you get exactly r 1's, output p to a file.
4. repeat from a) until you have a very large number of p's in your file.
5. Look at the cumulative histogram of the p's and clip 2.5% off each end to give you the 95% confidence interval.

Here is a very simple C++ program to do steps 1) to 4)

#include <stdio.h>
#include <stdlib.h>
#include <iostream.h>
#include "abmath.h"  // Provides drandom(), gives a double between 0
// and 1

main( int argc, char * argv[])
{
  int n, r, rr;
  const int N = 10000;
  char * tailPtr;
  double p;

  n = strtol( argv[ 1], &tailPtr, 10);
  r = strtol( argv[ 2], &tailPtr, 10);

  cerr << "(n,r)=(" << n << ", " << r << ")\n";

  for( int i=0; i < N;)
    {
      p = drandom();

      rr = 0;
      for( int j=0; j<n; j++)
        if( drandom() < p)
          rr++;

      if( rr == r)
        {
          cout << p << '\n';
          i++;
        }
    }
}