Skip to content

Instantly share code, notes, and snippets.

@ryamada22
Created October 16, 2018 09:13
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save ryamada22/4ceb0b001ccfd3a216581b29bfa54aa7 to your computer and use it in GitHub Desktop.
# tab 3 genotype count
my.hwe <- function(tab){
p.exp <- (2*tab[1]+tab[2])/(2*sum(tab))
g.exp <- c(p.exp^2,2*p.exp*(1-p.exp),(1-p.exp)^2)
N <- sum(tab)
chisq <- sum((tab-N*g.exp)^2/(N*g.exp))
return(chisq)
}
n.iter <- 10^4
chisq.stock <- rep(0,n.iter)
for(i in 1:n.iter){
p <-runif(1) * 0.6 + 0.2
N <- sample(500:1000,1)
obs <- sample(0:2,N,replace=TRUE,prob=g)
tab <- table(obs)
chisq <- my.hwe(tab)
chisq.stock[i] <- chisq
}
hist(chisq.stock)
k <- 1
chisq.theory <- rchisq(n.iter,df=k)
plot(sort(chisq.stock),sort(chisq.theory))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment