simulate Hardy-Weinberg Equilibrium using Markov chain

hwe_markov_chain.R

hwe_markov <- function(theta,g) {

#Taken from Ch7 question 7.9 from 
# E.A. Suess and B.E. Trumbo 
#Introduction to Probability and Gibbs Sampling with R 2010


# simulate Hardy Weinberg equlibrium using Markov chains

#Let the genotypes
#aa = 1, Aa = 2, and AA = 3 be the states of a process. At step 1, a female
#is of genotype aa, so that X1 = 1. At step 2, she selects a mate at random
#and produces one or more daughters, of whom the eldest is of genotype X2.
#At step 3, this daughter selects a mate at random and produces an eldest
#daughter of genotype X3, and so on

#number of generations is g
# theta is minor allele frequency

x = numeric(g) 
x[1] = 1

#the transition matrix:

# (theta, 1-theta, 0,
#  theta/2, .5, (1-theta)/2,
#  0, theta, (1-theta) )

#Simulation
for (i in 2:g)
{

if ( x[i-1] == 1 ) 
   x[i] = sample(1:3, 1, prob=c(theta, 1-theta, 0) )


if ( x[i-1] == 2 ) 
   x[i] = sample(1:3, 1, prob=c( (theta/2), .5, (1-theta)/2)  )


if ( x[i-1] == 3 ) 
   x[i] = sample(1:3, 1, prob=c(0, theta, 1-theta))

}

s <- summary(as.factor(x))/g # Table of proportions
s
#the stationary distribution should be (theta^2, 2*theta(1-theta), (1-theta)^2)

}