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)
  }