[R] table problems
Jim Lemon
bitwrit at ozemail.com.au
Wed Jun 12 09:42:32 CEST 2002
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);
}
More information about the R-help
mailing list