--- a
+++ b/weighted.statsOLD.R
@@ -0,0 +1,97 @@
+weighted.stats =
+  function (x, w, c) 
+  {
+    n = ncol(x) #number of samples
+    p = nrow(x) #number of genes
+    k = length(unique(c)) #number of class
+    c = as.integer(c)
+    WM = WS = matrix(0, p, k)
+    rownames(WS) = rownames(x)
+    colnames(WS) = unique(c)
+    c0 <- as.integer(min(c, na.rm = TRUE) - 1)
+    c <- as.integer(c) - c0
+    mk = NULL
+    
+    w.mean00 =
+      function (x, w) 
+      {
+        wm = NULL
+        
+        for (i in 1:p)
+        {
+          wm0 = sum(w[i,]*x[i,]) / sum(w[i,])
+          wm = c(wm, wm0)
+        }
+        return(wm)
+      }
+    
+    w.mean =
+      function (x, w, c) 
+      {
+        for (j in 1:k)
+        {
+          WM[,j] = w.mean00(x[,c == j], w[,c == j])
+        }
+        return(WM)
+      }
+    
+    w.sd =
+      function (x, w, c) 
+      {
+        w.sd00 =
+          function (x, w) 
+          {
+            ws = NULL
+            
+            w.sd0 =
+              function (x, w)
+              {
+                sumw = sum(w)
+                sumw.sq = sum(w)^2
+                w.sq = sum(w^2)
+                denom = sum(w * ((x - mean(x))^2))
+                sqrt((sumw * denom) / (sumw.sq - w.sq))
+              }
+            
+            for (i in 1:p)
+            {
+              ws0 = w.sd0(x[i,], w[i,])
+              ws = c(ws, ws0)
+            }
+            
+            return(ws)
+          }
+        
+        for (j in 1:k)
+        {
+          WS[,j] = w.sd00(x[,c == j], w[,c == j])
+        }
+        return(WS)
+      }
+    
+    WMEAN = w.mean00(x, w) #Overall weighted mean
+    delta = WMEAN.G = w.mean(x, w, c) #Weighted means for each group 
+    WSD.G = w.sd(x, w, c) #Weighted standard deviations for each group
+    WSD.POOLED = WSD.G
+    
+    for (i in 1:k)
+    {
+      WSD.POOLED[,i] = (table(c)[i]-1) * (WSD.POOLED[,i]^2)
+      mk[i] = sqrt((1 / table(c)[i]) + (1 / n))
+    }
+    
+    WSD.POOLED = sqrt(rowSums(as.data.frame(WSD.POOLED)) / (n - k))
+    s0 = median(WSD.POOLED)
+    
+    for (i in 1:k)
+    {
+      delta[,i] = (delta[,i] - WMEAN) / (mk[i]*(WSD.POOLED + s0))
+    }
+    
+    rownames(WMEAN) = rownames(WMEAN.G) = rownames(delta) = rownames(WSD.G) = rownames(WSD.POOLED) = rownames(x)
+    colnames(WMEAN.G) = colnames(delta) = colnames(WSD.G) = unique(c)
+    
+    stats = list(n = n, p = p, k = k, mk = mk, s0 = s0, delta = delta, WMEAN = WMEAN, WMEAN.G = WMEAN.G, WSD.G = WSD.G, WSD.POOLED = WSD.POOLED)
+    
+    return(stats)
+  }
\ No newline at end of file