## https://cmustatistics.github.io/data-repository/biology/dry-beans.html

rm(list=ls())


Multinomial = function(X,Y){
     ### multinomial logistic regression
     ### assume y takes values in {0,1,2,...,k}

     X = as.matrix(X)
     d = ncol(X)
     n = length(Y)
     y = unique(Y)
     k = length(y)

     beta = matrix(0,k-1,d+1)
     

     I = (1:n)[Y == y[1]]
     X0 = as.matrix(X[I,])
     n0 = nrow(X0)

     for(j in 1:(k-1)){
          J = (1:n)[Y == y[j+1]]
          X1 = as.matrix(X[J,])
          n1 = nrow(X1)

          XX = rbind(X0,X1)
          Z  = rep(1,n0+n1)
          Z[1:n0] = 0

          tmp = glm(Z ~ XX,family="binomial")$coef
 
          beta[j,] = tmp
          }


    one = rep(1,n)

    XX = cbind(one,X)

    Xbeta = matrix(0,n,k-1)

    for(j in 1:(k-1)){
        Xbeta[,j] = XX %*% as.matrix(beta[j,])
        }


    P = matrix(0,n,k)
    Pred = rep(0,n)

    for(i in 1:n){

        e = exp(Xbeta[i,])
        P[i,1] = 1/(1+sum(e))
        P[i,-1] = P[i,1]*e

        Pred[i] = y[which.max(P[i,])]
        }


     Confusion = table(factor(Y,levels=y),factor(Pred,levels=y))

     return(list(beta=beta,P=P,Pred=Pred,Confusion=Confusion))
     }



D = read.csv("dry-beans.csv",header=TRUE)
attach(D)

##  names(D)
##  [1] "Area"            "Perimeter"       "MajorAxisLength" "MinorAxisLength"
##  [5] "AspectRation"    "Eccentricity"    "ConvexArea"      "EquivDiameter"  
##  [9] "Extent"          "Solidity"        "roundness"       "Compactness"    
## [13] "ShapeFactor1"    "ShapeFactor2"    "ShapeFactor3"    "ShapeFactor4"   
## [17] "Class"

n = nrow(D)


Y = D$Class
X = D[,1:16]     

out = Multinomial(X,Y)
names(out)
## [1] "beta"      "P"         "Pred"      "Confusion" "tmp"

Confusion = out$Confusion

## Confusion
##            SEKER BARBUNYA BOMBAY CALI HOROZ SIRA DERMASON
##   SEKER     1924        5      0    0     0   60       38
##   BARBUNYA   863       93      0  178   161   26        1
##   BOMBAY     521        0      0    1     0    0        0
##   CALI      1588        1      0    6    31    4        0
##   HOROZ     1880        0      0    0    45    0        3
##   SIRA       184        2      0    0  1259  958      233
##   DERMASON   366        0      0    0   782   45     2353