Hi there,
I am trying to write a function to perform GOF test of data to a
zero-truncated Poisson distribution. I am facing 2 problems.
1) How can I obtain a frequency table for all values within the range of
observed values?
For instance if the observations are
obs <- c("A", "A", "A", "A", "B", "C", "C", "D", "E", "E", "F", "G",
"G", "H", "H", "H", "H")
## counts how many times indiviuals are counted 1, 2, 3, 4 times and so
(Fi <- table(table(obs)))
Fi
1 2 4
3 3 2
meaning that 3 individuals have been counted once, 3 counted twice and 2
counted 4 times. Applying table() here does not return the frequency of
individuals counted 3 times (which is 0 in this case). How can I achieve
such a contingency table where I could end up with something like
Fi
1 2 3 4
3 3 0 2
2) function to estimate lambda with MLE
trunpoismle <- function(xbar, tol = .Machine$double.eps^0.25, nmax = 10000)
{
lambda <- xbar
i <- 1
diff <- 1
while (abs(diff) > tol && i < nmax)
{
diff <- (lambda - xbar * (1 - exp(-lambda))) /
(1 - xbar * exp(-lambda))
lambda <- lambda - diff
i <- i + 1
}
if (i == nmax)
{
print("Iteration did not converge")
}
return(lambda)
}
3) Once I get there, I would like to apply a chi-square with (k-2) df, k
being the number of classes. Here the difficulty is to reduce the number
of classes so that a chi-square test may be applied. For instance:
Fi <- c(8, 5, 0, 4, 0, 1, 2, 0, 0, 1)
muHat <- trunpoismle(mean(Fi))
class <- 1 : 10
Exp <- muHat^class * exp(-muHat) / factorial(class)
Exp <- Exp[-1] / (1- Exp[1]) * sum(Fi)
# individuals seen 1, 2, 3, 4, 5, 6, 7, 8, 9, and 10 times
# This distribution maybe reduced to:
FiPooled <- c(8, 5, 4, 4)
# individuals seen 1, 2, 3-4 and > 4 times to apply chi-square tests
I would like to collapse the expected values the same way. Does anyone
have some tips or point out on a function I could have missed to pool
classes of a contingency table (or even to perfom the GOF test itself).
Thanks,
Roger
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
