[R] table problems

[Jim Lemon]

Robin Hankin wrote:
> 
>... I've been advised that the correct distribution for the number of
> sleeps in the 50 trees would be multinomial with n=50 and the
> p_i=12/50 for i=1 to 50.  I'll check out Feller when I get home
> tonight, but I suspect that M-B is the limiting case?
> 
I grabbed my thesis and refreshed my memory, and I think it will be what
you want. Here is a table for the discrete case with 50 cells and 12
trials.

N choices     probability     cumulative p
      1        2.048e-19       2.048e-19
      2        2.05421e-14     2.05423e-14
      3        4.16787e-11     4.16992e-11
      4        1.3844e-08      1.38857e-08
      5        1.43652e-06     1.45041e-06
      6        6.20311e-05     6.34815e-05
      7        0.00129369      0.00135717
      8        0.0141003       0.0154574
      9        0.0829514       0.0984088
     10        0.260324        0.358733
     11        0.403082        0.761815
     12        0.238185        1
 
Expected value = 10.76

What this is telling you is that if your possum is going back to the
same tree (which the ringtail at my place does) he will have to restrict
his choices to 8 or less to avoid a Type I error. Of course, your design
probably includes more than 1 possum, and you would have to sort out the
alpha for your number of subjects (and repeats?). I've appended the
source code for generating tables such as this.

Jim
-------------- next part --------------
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>

double Factorial(double n) {
 if(n > 1) return(n * Factorial(n-1));
 return(1);
}

double NChooseR(double n,double r) {
 double nfact;
 double rfact;
 double nminusrfact;

 nfact = Factorial(n);
 rfact = Factorial(r);
 nminusrfact = Factorial(n-r);
 return(nfact/(rfact * nminusrfact));
}

double PFilled(double n,double r,int f) {
 double psum=0;
 double combfact;
 int v;

 combfact = NChooseR(n,n-f);
 for(v=0;v<=f;v++)
  psum += pow(-1,v) * NChooseR(f,v) * pow((f-v)/n,r);
 return(combfact*psum);
}

double MBdist(double ncells,double nballs) {
 double psum;
 double ptotal = 0;
 double ev = 0;
 int filled = 1;

 while(filled <= nballs) {
  psum = PFilled(ncells,nballs,filled);
  ptotal += psum;
  fprintf(stdout,"     %2d        %-16g%g\n",filled,psum,ptotal);
  ev += psum * filled++;
 }
 if(ptotal < 0.99 || ptotal > 1.01)
  fprintf(stdout,"Warning: total p is %f\n",ptotal);
 return(ev);
}
 
int main(int argc,char *argv[]) {
 if(argc == 3) {
  fprintf(stdout,"N choices     probability     cumulative p\n");
  printf("\nExpected value = %.2f\n",MBdist(atof(argv[1]),atof(argv[2])));
  return(0);
 }
 return(1);
}