[R] Getting all possible contingency tables

(Ted Harding) Ted.Harding at wlandres.net
Sun Dec 2 16:02:48 CET 2012


On 02-Dec-2012 14:17:20 Christofer Bogaso wrote:
> Thanks Ted (and other) for your suggestion. Here I have implemented 
> following:
> 
> Tab <- matrix(c(8, 10, 12, 6), nr = 2)
> 
> Simu_Number <- 50000
> Tab_Simulate <- vector("list", length = Simu_Number)
> for (i in 1:Simu_Number) {
>          Tab_Simulate[[i]] <- matrix(rmultinom(1, sum(Tab), rep(0.25, 
> 4)), nrow = 2)   ## All Cells have equal probability
>      }
> Sample_ChiSq <- sapply(Tab_Simulate, function(x) {
>                                  Statistic <-
> sum((chisq.test(as.table(x))$observed - 
> chisq.test(as.table(x))$expected)^2/chisq.test(as.table(x))$expected)
>                                  return(Statistic)
>                              })
> 
> length(Sample_ChiSq[Sample_ChiSq < qchisq(0.95, 1)])/Simu_Number
> 
> hist(Sample_ChiSq, freq = FALSE)
> lines(dchisq(seq(min(Sample_ChiSq), max(Sample_ChiSq), by = 0.5), 1))
> 
> 
> However I think I am making some serious mistake as histogram did not 
> match the density curve.
> 
> Can somebody help me where I am making mistake?
> 
> Thanks and regards,
> [the remainder (copies of previous posts) snipped]

The reasons for the mis-match are:

A: that you have put the curve in the wrong place, by not
supplying x-coordinates to lines(), so that it plots its
points at x = 1,2,3,4,...

B: that you need to multiply the plotted density by the width
of the histogram cells, so as to match the density curve to the
discrete density of the histogram. It will also then look better
when the chis-squared curve is plotted at the mid-points of the cells.

Hence, try something like:

hist(Sample_ChiSq, freq = FALSE, breaks=0.5*(0:40))
x0 <- 0.25+0.5*(0:20)
lines(x0,dchisq(x0,1))

Hoping this helps,
Ted.

-------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at wlandres.net>
Date: 02-Dec-2012  Time: 15:02:45
This message was sent by XFMail




More information about the R-help mailing list