# R coding of genomic selection from Meuwissen et al. (2001)

# Parameters
nmarkers <- 10      #number of markers
nrecords <- 325     #number of records
numit    <- 1000   #number of iterations

# Read in data
x <- matrix(scan("c:/texasam2011/practicals/day5/realDataExample/xvec_day4.inp"),ncol=nmarkers,byrow=TRUE)
y <- matrix(scan("c:/texasam2011/practicals/day5/realDataExample/yvec_day4.inp"),byrow=TRUE)

# Storage vectors and matricies
gStore <- array(0,c(numit,nmarkers))
gvarStore <- array(0,c(numit,nmarkers))
vareStore <- array(0,c(numit))
muStore <- array(0,c(numit))
ittstore <- array(0,c(numit))

# Algorithm
# Step 1.  Initialization
g <- array(0.01,c(nmarkers)) 

mu <- 0.1
gvar <- array(0.1,c(nmarkers))
ones <- array(1,c(nrecords))
e <- array(0,c(nrecords))

for (l in 1:numit) {

# Step 2.  Sample the gvar from the inverse chi square posterior
      for (j in 1:nmarkers) {
# Draw from inverse chi squared distribution
#        gvar[j] <- (0.002+g[j]*g[j])/rchisq(1,4.012+1)  # Meuwissen et al. (2001) prior
         gvar[j] <- (g[j]*g[j])/rchisq(1,1+1)        # Xu (2003) prior
#        gvar[j] <- (g[j]*g[j])/rchisq(1,0.998)    # Te Braak et al. (2006) prior
      }

# Step 3.  Sample vare from an inverse chi-square posterior   
      e <- y - x%*%g - mu 
      vare <- (t(e)%*%e/(nrecords-2))/rchisq(1,nrecords-2)
            
# Step 4 Sample the mean from a normal posterior 
      mu <- rnorm(1,(t(ones)%*%y - t(ones)%*%x%*%g)/nrecords,sqrt(vare/nrecords))

# Step5 Sample the g from a normal distribution
      z <- array(0,c(nrecords))
      gold <- g
      for (j in 1:nmarkers) {
       gtemp <- g
       gtemp[j] <- 0
       for (i in 1:nrecords) {
        z[i] <- x[i,j]
       }
       mean <- ( t(z)%*%y-t(z)%*%x%*%gtemp-t(z)%*%ones*mu ) / (t(z)%*%z+vare/gvar[j])
       g[j] <- rnorm(1,mean,sqrt(vare/(t(z)%*%z+vare/gvar[j])))
      }

# Store values for each itteration
      for (j in 1:nmarkers) {
        gStore[l,j] <- g[j]
        gvarStore[l,j] <- gvar[j]
      }
      vareStore[l] <- vare
      muStore[l] <- mu 
      ittstore[l] <- l
}
meanSNPeffect <- array(0,c(nmarkers))
for (j in 1:nmarkers) {
 meanSNPeffect[j] <- mean(gStore[100:1000,j])
}
meanSNPeffect
 


par(mfcol=c(5,4))
par(mai=c(0.1,0.1,0.1,0.5)+0.2)

for (i in 1:nmarkers) {
 plot(ittstore[1:1000],gStore[1:1000,i])
}
for (i in 1:nmarkers) {
 plot(density(gStore[100:1000,i]))
}