cdg-snl.R

cdg <- read.csv("cdgTable.csv2", header = T, na.strings = T)
gene1 <- read.csv("gene-1.snp.csv", header = T, row.names = 1)

snps <- t(gene1)
cdg.mean <- tapply(cdg$logcdg, cdg$strains, mean)
t <- as.matrix(ifelse(cdg.mean[row.names(snps)] < -1, 0,1 ))
snps <- cbind(t, snps)

w <- runif(32, 1e-3, 1e-2)
b <- runif(1)

eta <- 0.1

weights <- matrix(0,nrow = 1000, ncol = 32)
accuracy <- numeric()

for(epoch in 1:1000){
  a <- snps[,2:33] %*% w
  y <- 1/(1+exp(-a-b))
  e <- snps[,1] - y
  w <- w - eta * -colSums(snps[,2:33] * e[,1])
  b <- b - eta * sum(-e)
  for(k in 1:ncol(weights)){
   weights[epoch,k] <- w[k]
  }
  accuracy <- c(accuracy, length(which(round(y) == t))/30)
}

plot(accuracy, type = "p")

plot(NA, xlim = c(1,1000), ylim = c(-8.78753,18.08474))
for(i in 1:32){
  lines(x = 1:1000, y = weights[,i], col=sample(rainbow(150),1), lty=1)
}