--- a
+++ b/classnb.R
@@ -0,0 +1,43 @@
+classnb <- function (x, y, xte = NULL, beta = 1,  type = c("mle", 
+    "deseq", "quantile"), prior = NULL) 
+{
+
+    if (is.null(prior)) 
+        prior <- rep(1/length(unique(y)), length(unique(y)))
+
+    null.out <- NullModel(x, type = type)
+    ns <- null.out$n
+    nste <- NullModelTest(null.out, x, xte, type = type)$nste
+
+    uniq <- sort(unique(y))
+
+        ds <- GetDnb(ns, x, y, beta)
+
+        #disperhat=rep(truephi,ncol(nste))
+        phihat <- as.numeric(disperhat)
+        discriminant <- matrix(NA, nrow = nrow(xte), ncol = length(uniq))
+
+        for (k in 1:length(uniq)) {
+
+            for(l in 1:nrow(xte))   {
+
+                 dstar = ds[k,]
+                 part2=1+nste[l,]*dstar*phihat 
+                 part1=dstar/part2 
+
+                 discriminant[l, k]<- sum(xte[l,]*log(part1))-sum((1/phihat)*log(part2))+log(prior[k])
+
+             }
+         }
+        save <- list(ns = ns, nste = nste, ds = ds, discriminant = discriminant, 
+            ytehat = uniq[apply(discriminant, 1, which.max)], 
+            x = x, y = y, xte = xte, 
+            type = type)
+
+        return(save)
+
+}
+
+
+
+