Example R programs and commands 23. Poisson density and tests for randomness # All lines preceded by the "#" character are my comments. # All other left-justified lines are my input. # All other indented lines are the R program output. # POISSON probability density # # Example: probability of x=5 counts when the mean is lambda=7. x<-5 # event: count lambda<-7 # situation: expected count dpois(x, lambda) # probability of count x when the mean is lambda # Complete list of probabilities for x=0,1,2,...,(big): dpois(0:20, lambda) # probabilities of counts x=0,1,...,20 # Cumulative Poisson density function: ppois(x, lambda) # probability of x or fewer counts ppois(x, lambda, lower.tail=FALSE) # probability of more than x counts 1-ppois(x, lambda) # ...same result as the previous line # POISSON randomness tests # # EG: test if the following frequency table has a Poisson density. # Raisin count: 0 1 2 3 4 5 6 or more # Frequency: 5 12 18 14 5 1 0 # rcnt<-0:6 # Table contains frequencies of counts 0,1,...,6 freq<-scan() 5 12 18 14 5 1 0 lasti<-length(rcnt) # last index in the rcnt[] and freq[] arrays. rcbar<- sum(rcnt*freq)/sum(freq); rcbar # average raisin count per slice phat<-dpois(rcnt,lambda=rcbar) # Poisson density for the observed raisin counts # ... piling the probabilities for all counts bigger than # the maximum count with a nonzero frequency into the last value: phat[lasti] <- ppois(rcnt[lasti-1], lambda=rcbar, lower.tail=FALSE) # NOTE: the expression finds the actual last count with nonzero frequency: max(rcnt*(freq>0)) # Perform the goodness-of-fit test after correcting the number of # degrees of freedom to be lasti-2: pchisq(chisq.test(freq,p=phat)\$statistic,df=lasti-2,lower.tail=FALSE) # OUTPUT: # X-squared # 0.4656001 # Warning message: # Chi-squared approximation may be incorrect in: chisq.test(freq, p = phat) # # The p-value is 0.4656001. # Conclusion: DO NOT REJECT H0: raisins are randomly distributed by slice. # Example with decidedly NON-RANDOM raisin distribution. # Raisin count: 0 1 2 3 4 5 6 or more # Frequency: 0 0 0 0 0 50 0 # rcnt<-0:6 # Table contains frequencies of counts 0,1,...,6 freq<-scan() 0 0 0 0 0 50 0 lasti<-length(rcnt) # last index in the rcnt[] and freq[] arrays. rcbar<- sum(rcnt*freq)/sum(freq); rcbar # average raisin count per slice phat<-dpois(rcnt,lambda=rcbar) # Poisson density for the observed raisin counts # ... piling the probabilities for all counts bigger than # the maximum count with a nonzero frequency into the last value: phat[lasti] <- ppois(rcnt[lasti-1], lambda=rcbar, lower.tail=FALSE) # Perform the goodness-of-fit test after correcting the number of # degrees of freedom to be lasti-2: pchisq(chisq.test(freq,p=phat)\$statistic,df=lasti-2,lower.tail=FALSE) # OUTPUT: # X-squared # 9.27622e-49 # Warning message: # Chi-squared approximation may be incorrect in: chisq.test(freq, p = phat) # # The p-value is 9.27622e-49 # Conclusion: REJECT H0 for HA: raisins are NOT randomly distributed by slice.