.packageName <- "multinomRob"
#
#  multinomRob
#
#  Walter R. Mebane, Jr.
#  Cornell University
#  http://macht.arts.cornell.edu/wrm1/
#  wrm1@macht.arts.cornell.edu
#
#  Jasjeet Singh Sekhon 
#  Harvard University
#  http://jsekhon.fas.harvard.edu/
#  jsekhon@fas.harvard.edu
#
#  $Id: datamanip.R,v 1.5 2004/02/15 23:35:30 wrm1 Exp $
#

#Mapping from xvec to the beta.vector and back again
#forward==TRUE from xvec TO beta.vector
#forward==FALSE from beta.vector TO xvec
mnl.xvec.mapping <- function (forward=TRUE,base.xvec,work.xvec,beta.vector,
                             ncats,nvars.total) 
{
  n.ones <- sum(base.xvec==1);  # indicates unique parameters
  n.mults <- sum(base.xvec>1);  # indicates parameters constrained equal
  if (forward) {
    p <- 0;
    if (n.ones>0) {
      for (j in 1:ncats) {
        for (i in 1:nvars.total) {
          if (base.xvec[i,j] == 1) {
            p <- p + 1;
            beta.vector[p] <- work.xvec[i,j];
          } #end of if
        } #end of i
      } #end of j
    }
    if (n.mults>0) {
      nidxvals <- length(idxvals <- sort(unique(base.xvec[base.xvec>1])));
      for (k in 1:nidxvals) {
        for (j in 1:ncats) {
          if (any(jktest <- base.xvec[,j]==idxvals[k])) {
            beta.vector[p+k] <- work.xvec[jktest,j][1];
            break;
          }
        }
      }
    }
    return(beta.vector);
  } else {
    p <- 0;
    if (n.ones>0) {
      for (j in 1:ncats) {
        for (i in 1:nvars.total) {
          if (base.xvec[i,j] == 1) {
            p <- p + 1;
            work.xvec[i,j] <- beta.vector[p];
          } #end of if
        } #end of i
      } #end of j
    }
    if (n.mults>0) {
      nidxvals <- length(idxvals <- sort(unique(base.xvec[base.xvec>1])));
      for (k in 1:nidxvals) {
        for (j in 1:ncats) {
          if (any(jktest <- base.xvec[,j]==idxvals[k])) {
            work.xvec[jktest,j] <- beta.vector[p+k];
          }
        }
      }
    }
    return(work.xvec);
  } # end of else
} #end of mnl.xvec.mapping


############################################################################
## Create the jacstack (from tanh)                                     #
############################################################################    

# jacstack.function:  arrange regressors for computing the Jacobian matrix
##   produces jacstack:  array of regressors,
##      dim(jacstack) = c(n observations, n UNIQUE parameters, n categories)
jacstack.function <- function(X,tvars.unique,xvec) {
  xdim  <- dim(X)
  obs   <- xdim[1]
  nvars <- xdim[2]
  ncats <- xdim[3]

  jacstack <- array(0,dim=c(obs,tvars.unique,ncats));
  n.ones <- sum(xvec==1);  # indicates unique parameters
  n.mults <- sum(xvec>1);  # indicates parameters constrained equal
  itu <- 0;
  if (n.ones>0) {
    for (j in 1:ncats) {
      nxj <- sum(xvec[,j]==1);
      if (nxj>0) {
        jacstack[1:obs, itu + 1:nxj, j] <- X[rep(TRUE,obs),xvec[,j]==1,j];
      }
      itu <- itu + nxj;
    }
  }
  if (n.mults>0) {
    nxvecrows <- dim(xvec)[1];
    nidxvals <- length(idxvals <- sort(unique(xvec[xvec>1])));
    for (k in 1:nidxvals) {
      for (j in 1:ncats) {
        if (any(jktest <- xvec[,j]==idxvals[k])) {
          kidx <- (1:nxvecrows)[jktest];  # xvec[,j] rows matching constraint k
          nxj <- sum(jktest);
          for (kk in 1:nxj) {
            jacstack[1:obs, itu + k, j] <-
              jacstack[1:obs, itu + k, j] + X[rep(TRUE,obs),kidx[kk],j];
          }
        }
      }
    }
  }
  return(jacstack)
} #jacstack.function

# jacstack.singles:  check jacstack array for regressors with a distinct value
#                    at only one observation
jacstack.singles <- function(jacstack) {
  nunique <- dim(jacstack)[2];

  jsingle <- rep(FALSE, nunique);
  for (i in 1:nunique) {
    if (length(unique(c(jacstack[,i,]))) == 2) {
      jsingle[i] <- ifelse(any(table(c(jacstack[,i,])) == 1), TRUE, FALSE)
    }
  }
  return(jsingle)
} #jacstack.singles.function
#
#  multinomRob
#
#  Walter R. Mebane, Jr.
#  Cornell University
#  http://macht.arts.cornell.edu/wrm1/
#  wrm1@macht.arts.cornell.edu
#
#  Jasjeet Singh Sekhon 
#  Harvard University
#  http://jsekhon.fas.harvard.edu/
#  jsekhon@fas.harvard.edu
#
#  $Id: genoudRob.R,v 1.4 2004/03/04 02:08:17 wrm1 Exp $
#
###################################
#New Front End for Genoud, with tuned defaults
###################################

#sets genoud.parms defaults
genoudParms  <- function(genoud.parms)
  {
    #set user controlled defaults
    if (is.null(genoud.parms$pop.size))
      genoud.parms$pop.size  <- 1000;

    if (is.null(genoud.parms$max.generations))
      genoud.parms$max.generations  <- 100;
    
    if (is.null(genoud.parms$wait.generations))
      genoud.parms$wait.generations  <- 10;
    
    if (is.null(genoud.parms$hard.generation.limit))
      genoud.parms$hard.generation.limit  <- FALSE;

    #this is redundant, but maintains clarity  
    if (is.null(genoud.parms$MemoryMatrix))
      genoud.parms$MemoryMatrix  <- NULL;
  
    if (is.null(genoud.parms$Debug))
      genoud.parms$Debug  <- FALSE ;

    #this is redundant, but maintains clarity   
    if (is.null(genoud.parms$Domains))
      genoud.parms$Domains  <- NULL;

    if (is.null(genoud.parms$scale.domains))
      genoud.parms$scale.domains  <- 10;
    
    if (is.null(genoud.parms$boundary.enforcement))
      genoud.parms$boundary.enforcement  <- 0;
    
    if (is.null(genoud.parms$solution.tolerance))
      genoud.parms$solution.tolerance  <- 0.0000001;
    
    if (is.null(genoud.parms$BFGS))
      genoud.parms$BFGS  <- TRUE;
  
    if (is.null(genoud.parms$unif.seed))
      genoud.parms$unif.seed  <- 812821;
    
    if (is.null(genoud.parms$int.seed))
      genoud.parms$int.seed  <- 53058;
    
    if (is.null(genoud.parms$print.level))
      genoud.parms$print.level  <- 0;
    
    if (is.null(genoud.parms$share.type))
      genoud.parms$share.type  <- 0;
    
    if (is.null(genoud.parms$instance.number))
      genoud.parms$instance.number  <- 0;
    
    if (is.null(genoud.parms$output.path))
      genoud.parms$output.path  <- "stdout";
    
    if (is.null(genoud.parms$output.append))
      genoud.parms$output.append  <- FALSE;
    
    if (is.null(genoud.parms$project.path))
      genoud.parms$project.path  <- "/dev/null";
    
    if (is.null(genoud.parms$P1))
      genoud.parms$P1  <- 50;
    
    if (is.null(genoud.parms$P2))
      genoud.parms$P2  <- 50;
    
    if (is.null(genoud.parms$P3))
      genoud.parms$P3  <- 50;
    
    if (is.null(genoud.parms$P4))
      genoud.parms$P4  <- 50;
    
    if (is.null(genoud.parms$P5))
      genoud.parms$P5  <- 50;
    
    if (is.null(genoud.parms$P6))
      genoud.parms$P6  <- 50;
    
    if (is.null(genoud.parms$P7))
      genoud.parms$P7  <- 50;
    
    if (is.null(genoud.parms$P8))
      genoud.parms$P8  <- 50;
    
    if (is.null(genoud.parms$P9))
      genoud.parms$P9  <- 0  ;

    return(genoud.parms);
  } #end genoudParms


genoudRob <- function(fn,nvars,starting.values,genoud.parms)
{
  #set static defaults
  max  <- FALSE
  gradient.check  <- FALSE
  data.type.int  <- FALSE
  hessian  <- FALSE  
  roptim <- TRUE;

  #load up genoud.parms
  pop.size  <- genoud.parms$pop.size;
  max.generations  <- genoud.parms$max.generations;
  wait.generations  <- genoud.parms$wait.generations;
  hard.generation.limit  <- genoud.parms$hard.generation.limit;
  MemoryMatrix  <- genoud.parms$MemoryMatrix;
  Debug  <- genoud.parms$Debug;
  Domains  <- genoud.parms$Domains;
  scale.domains  <- genoud.parms$scale.domains;
  boundary.enforcement  <- genoud.parms$boundary.enforcement;
  solution.tolerance  <- genoud.parms$solution.tolerance;
  BFGS  <- genoud.parms$BFGS;
  unif.seed  <- genoud.parms$unif.seed;
  int.seed  <- genoud.parms$int.seed;
  print.level  <- genoud.parms$print.level;
  share.type  <- genoud.parms$share.type;
  instance.number  <- genoud.parms$instance.number;
  output.path  <- genoud.parms$output.path;
  output.append  <- genoud.parms$output.append;
  project.path  <- genoud.parms$project.path;
  P1  <- genoud.parms$P1;
  P2  <- genoud.parms$P2;
  P3  <- genoud.parms$P3;
  P4  <- genoud.parms$P4;
  P5  <- genoud.parms$P5;
  P6  <- genoud.parms$P6;
  P7  <- genoud.parms$P7;
  P8  <- genoud.parms$P8;
  P9  <- genoud.parms$P9;

  #we always have starting, but leave this check in. 
  #do we have starting values?
  if (is.null(starting.values)) {
    nStartingValues <- 0;
    parm.vec  <- rep(1,nvars)
  }
  else {
    nStartingValues <- 1;
    parm.vec  <- starting.values;
  }

  # let's create the Domains if none have been passed.
  if (!(is.matrix(Domains)))
    {
      Domains <- matrix(nrow=nvars, ncol=2);
      for (i in 1:nvars)
        {
          Domains[i,1] <- parm.vec[i] - abs(parm.vec[i])*scale.domains;
          Domains[i,2] <- parm.vec[i] + abs(parm.vec[i])*scale.domains;
        } # end of for loop
    } # end of Domains if
  
  #MemoryMatrix
  if (is.null(MemoryMatrix)) {
    MemoryMatrix <- TRUE;

    if (nvars > 20) {
      if (print.level > 0) {
        cat("\nWARNING: Since the number of parameters is greater than 20,\nWARNING: MemoryMatrix has been turned off by default.\nWARNING: You may turn it on using the MemoryMatrix flag.\nWARNING: This option increases speed at the cost of extra memory usage.\n\n")
        MemoryMatrix <- FALSE;
      }
    }
  }

  #set output.type
  if (output.path=="stdout")
    {
      output.type <- 0;
    }
  else
    {
      if (output.append)
        {
          output.type <- 2;
        }
      else
        {
          output.type <- 1;
        }
    }

  # create the P vector
  P <- vector(length=9, mode="numeric");
  P[1] <- P1; P[2] <- P2; P[3] <- P3; P[4] <- P4;
  P[5] <- P5; P[6] <- P6; P[7] <- P7; P[8] <- P8;
  P[9] <- P9;

  # has the user provided any seeds?
  if (unif.seed==812821 && int.seed==53058)
    provide.seeds <- FALSE
  else
    provide.seeds <- TRUE;

  if (max==FALSE)
        {
          g.scale <- 1;
        }
  else
    {
      g.scale <- -1;
    }

  #optim st
  genoud.optim.wrapper101 <- function(foo.vals)
    {
      ret <- optim(foo.vals, fn=as.function(fn), method="BFGS",
                  control=list(fnscale=g.scale));
      return(c(ret$value,ret$par));
    } # end of genoud.optim.wrapper101


  gout <- .Call("rgenoud", as.function(fn), new.env(),
                as.integer(nvars), as.integer(pop.size), as.integer(max.generations),
                as.integer(wait.generations),
                as.integer(nStartingValues), as.vector(starting.values),
                as.vector(P), as.matrix(Domains),
                as.integer(max), as.integer(gradient.check), as.integer(boundary.enforcement),
                as.double(solution.tolerance), as.integer(BFGS), as.integer(data.type.int),
                as.integer(provide.seeds), as.integer(unif.seed), as.integer(int.seed),
                as.integer(print.level), as.integer(share.type), as.integer(instance.number),
                as.integer(MemoryMatrix), as.integer(Debug),
                as.character(output.path), as.integer(output.type), as.character(project.path),
                as.integer(hard.generation.limit),
                as.function(genoud.optim.wrapper101), as.integer(roptim),
                PACKAGE="rgenoud");

  if (hessian==TRUE)
    {
      hes <- optim(gout[5:(nvars+4)], fn, method="BFGS", hessian=TRUE,
                  control=list(fnscale=g.scale));
      
      hes <- hes$hessian;

      ret <- list(value=gout[1], generations=gout[2], peakgeneration=gout[3], popsize=gout[4],
                 par=gout[5:(nvars+4)], gradients=gout[(nvars+5):(nvars+nvars+4)],
                 operators=gout[(nvars+nvars+5):(nvars+nvars+9+4)],
                 hessian=hes);
    }
  else
    {
      ret <- list(value=gout[1], generations=gout[2], peakgeneration=gout[3], popsize=gout[4],
                 par=gout[5:(nvars+4)], gradients=gout[(nvars+5):(nvars+nvars+4)],
                 operators=gout[(nvars+nvars+5):(nvars+nvars+9+4)]);
    }

  return(ret);
} #end of genoudRob()


    
#
#  multinomRob
#
#  Walter R. Mebane, Jr.
#  Cornell University
#  http://macht.arts.cornell.edu/wrm1/
#  wrm1@macht.arts.cornell.edu
#
#  Jasjeet Singh Sekhon 
#  Harvard University
#  http://jsekhon.fas.harvard.edu/
#  jsekhon@fas.harvard.edu
#
#  $Id: multinomLQD.R,v 1.2 2004/03/04 02:08:17 wrm1 Exp $
#

##
## Multinomial lqd functions.
##
## MORE IN C CODE

## fit.multinomial.C.lqd2:  compute the LQD interquartile difference spread estimate
##
fit.multinomial.C.lqd2 <- function(foo,X,Y,Ypos,xvec,tvec,ncats,nvars,nvars.unique,obs,TotalY) {
  tmp.vec     <-   mnl.xvec.mapping(forward=FALSE,
                                    xvec,
                                    tvec,
                                    foo,
                                    ncats,
                                    nvars);  
  
  if (all(Ypos)) {
    return(.Call("original_fit_lqd2",
                 as.integer(obs),
                 as.integer(ncats),
                 as.integer(nvars),
                 as.integer(nvars.unique),
                 as.real(tmp.vec),
                 as.real(Y),
                 as.real(X),
                 as.real(TotalY),
                 PACKAGE="multinomRob")
           );
  }
  else {
    lqd2 <- function(r,nparms) {
      obs <- length(r)
      h <- ceiling((obs+nparms)/2);
      hidx <- (h*(h-1))/2;
      dif <- abs(outer(r,r,"-")[outer(1:obs,1:obs,">")])
    #  qrt <- sort(dif)[hidx]  * 2.21914446599;
      qrt  <-  kth.smallest(dif, length(dif), hidx)* 2.21914446599;  
      return(qrt)
    }

    Y[!Ypos] <- 0;
    sres <- resfunc.lqd2(Y, Ypos, X, tmp.vec);
    return( lqd2(c(sres[!is.na(sres)]),nvars.unique) );
  }
} #end of fit.multinomial.C.lqd2

#fit.multinomial.lqd2 <- function(foo,X,Y,xvec,tvec,ncats,nvars,nvars.unique) 
#{
#  tmp.vec     <- mnl.xvec.mapping(forward=FALSE,xvec,tmp.vec,foo,
#                                  ncats,nvars);
#  y.prob      <- mnl.probfunc(Y,X,tmp.vec);
#
#  Sres.raw    <- res.std(Y,TotalY,y.prob)
#
#  fit.lqd2 <- lqd2(Sres.raw,nvars.unique);
#  return(fit.lqd2);
#} #end of fit.multinomial.lqd2

## resfunc.lqd2:  orthogonalized and standardized (for multinomial covariance) resids
resfunc.lqd2 <- function(Y, Ypos, Xarray, tvec) {
  if (all(Ypos)) {
    r <- res.std(Y, c(Y %*% rep(1,dim(Y)[2])), mnl.probfunc(Y, Ypos, Xarray, tvec));
  }
  else {
    nobs <- dim(Y)[1];
    ncats <- dim(Y)[2];
    r <- matrix(NA, nobs, ncats-1);
    phat <- mnl.probfunc(Y, Ypos, Xarray, tvec);
    hasall <- apply(Ypos, 1, sum) == ncats;
    nobsall <- sum(hasall);
    if (nobsall > 0) {
      Yuse <- matrix(Y[hasall,], nobsall, ncats);  # in case nobsall == 1
      puse <- matrix(phat[hasall,], nobsall, ncats);
      r[hasall,] <- res.std(Yuse, c(Yuse %*% rep(1,ncats)), puse);
    }
    hasless <- (1:nobs)[!hasall];
    for (i in hasless) {  # orthostd resids go into r[i,1:(nlesscats-1)]
      usecats <- Ypos[i,];
      nlesscats <- sum(usecats);
      Yuse <- matrix(Y[i,usecats], 1, nlesscats);
      puse <- matrix(phat[i,usecats], 1, nlesscats);
      r[i,1:(nlesscats-1)] <- res.std(Yuse, c(Yuse %*% rep(1,nlesscats)), puse);
    }
  }
  return( r );
}
#
#  multinomRob
#
#  Walter R. Mebane, Jr.
#  Cornell University
#  http://macht.arts.cornell.edu/wrm1/
#  wrm1@macht.arts.cornell.edu
#
#  Jasjeet Singh Sekhon 
#  Harvard University
#  http://jsekhon.fas.harvard.edu/
#  jsekhon@fas.harvard.edu
#
#  $Id: multinomMLE.R,v 1.8 2004/02/19 02:13:11 wrm1 Exp $
#
## multinomMLE:  maximum likelihood estimator for grouped multinomial GLM, with overdispersion
##   Y:  matrix of (overdispersed and contaminated) multinomial counts
##   Ypos:  matrix indicating which in Y are counts (TRUE) and which are not (FALSE).
##   Xarray:  array of regressors,
##      dim(Xarray) = c(n observations, n parameters, n categories)
##   xvec:  vector to indicate all the coefficient parameters in the model
##      (parms by ncats):
##      It has a 1 for an estimated parameter and a 0 otherwize.
##      example:
##      > xvec
##           [,1] [,2] [,3] [,4] [,5]
##      [1,]    1    1    1    1    0
##      [2,]    1    1    1    1    0
##      [3,]    1    1    1    1    0
##      [4,]    1    1    1    1    0
##   tvec: parms by ncats matrix (matrix with LQD estimates):
##      example:
##      > tvec
##                          Buchanan        Nader     Gore     Bush Other
##      int               -0.1641034    1.0735560 3.363641 4.151853     0
##      p(r,dg,d,r)96      2.4413780    0.3269827 3.207104 1.676676     0
##      p(cr,g,cd,cr)RV00 15.5333800 1149.4130000 2.039766 1.761392     0
##      pCuban            -6.4083750    0.3546630 1.966287 2.795598     0
##   jacstack:  array of regressors,
##      dim(jacstack) = c(n observations, n UNIQUE parameters, n categories)
multinomMLE <- function(Y, Ypos, Xarray, xvec, 
                        jacstack, itmax=100,
                        xvar.labels, choice.labels, print.level=0) {
  ## probfunc: matrix of estimated probabilities
  probfunc <-
    function(Y, Ypos, Xarray, tvec) {
      nobs <- dim(Y)[1]
      ncats <- dim(Y)[2]
      eta <- matrix(0,nobs,ncats)
      for (j in 1:ncats) {
        useobs <- Ypos[,j];
        eta[useobs,j] <- exp(Xarray[useobs,,j] %*% tvec[,j]);
      }
      return( c(1/(eta %*% rep(1,ncats))) * eta )
    }
  ## scorefunc:  score matrix
  scorefunc <-
    function(Ypos, nobs, nparms, N, presmat, jacstack) {
      scoremat <- matrix(0,nparms,nobs);
      for (i in 1:nobs) {
        usecats <- Ypos[i,];
        scoremat[,i] <- N[i] * presmat[i,usecats] %*% t(jacstack[i,,usecats]) ;  ## unweighted
      }
      return( scoremat )
    }
  ## hessianfunc:  hessian matrix
  hessianfunc <-
    function(Ypos, nobs, nparms, phat, N, jacstack) {
      H <- matrix(0,nparms,nparms)
      for (i in 1:nobs) {
        usecats <- Ypos[i,];
        pvec <- phat[i,usecats];
        wpvmat <- diag(pvec)-outer(pvec,pvec);  ## unweighted
        H0 <- N[i] * wpvmat;
        H <- H + jacstack[i,,usecats] %*% H0 %*% t(jacstack[i,,usecats])
      }
      return( H )
    }
  ## resfunc:  orthogonalized and standardized (for multinomial covariance) resids
  resfunc <-
    function(Y, Ypos, Xarray, tvec) {
      if (all(Ypos)) {
        r <- res.std(Y, c(Y %*% rep(1,dim(Y)[2])), probfunc(Y, Ypos, Xarray, tvec));
      }
      else {
        nobs <- dim(Y)[1];
        ncats <- dim(Y)[2];
        r <- matrix(0, nobs, ncats-1);
        phat <- probfunc(Y, Ypos, Xarray, tvec);
        hasall <- apply(Ypos, 1, sum) == ncats;
        nobsall <- sum(hasall);
        if (nobsall > 0) {
          Yuse <- matrix(Y[hasall,], nobsall, ncats);  # in case nobsall == 1
          puse <- matrix(phat[hasall,], nobsall, ncats);
          r[hasall,] <- res.std(Yuse, c(Yuse %*% rep(1,ncats)), puse);
        }
        hasless <- (1:nobs)[!hasall];
        for (i in hasless) {
          usecats <- Ypos[i,];
          nlesscats <- sum(usecats);
          ocats <- 1:(nlesscats-1);
          Yuse <- matrix(Y[i,usecats], 1, nlesscats);
          puse <- matrix(phat[i,usecats], 1, nlesscats);
          r[i,ocats] <- res.std(Yuse, c(Yuse %*% rep(1,nlesscats)), puse);
        }
      }
      return( r );
    }
  ## check convergence
  converged <-
    function(bnew,bold) {
      return( sqrt(sum((bnew-bold)^2)) < 1e-6*(sqrt(sum(bold^2)) + 1e-4) )
    }
  converged2 <-
    function(snew,sold) {
      return( sqrt(sum((snew-sold)^2)) <= 1e-8*sqrt(sum(sold^2)) )
    }

  ## begin data computations
  Y[!Ypos] <- 0;  # ensure noncounts are set to zero, for convenience
  ncats <- dim(Y)[2]
  tvec  <- xvec
  ## begin definition of variables used in GNstep that do not change over iterations
  nobs <- dim(Y)[1]
  ncats <- dim(Y)[2]
  catidx <- 1:(ncats-1);
  tvars.total <- dim(Xarray)[2]
  tvunique <- dim(jacstack)[2]
  mvec <- c(Y %*% rep(1,dim(Y)[2]));
  propmat <- Y / mvec;  ## transform observed counts to proportions
  ## end definition of variables used in GNstep that do not change over iterations

  #starting values
  bvec <- rep(1,tvunique)

  LogLik <-
    function(Y,Ypos,ipmatS,mvecS) {
      LLu <- -sum(Y[Ypos] * log(ipmatS[Ypos]));  ## negative loglikelihood
      ##   print(paste("unweighted:",LLu))
      return(list(LLu=LLu));
    }

  ## Newton algorithm given data
  GNstep <-
    function(bvec,sigma2GN,itmaxGN=100) {
      tvec <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,bvec, ncats,tvars.total);
      itersGN <- 0;
      ## patterned after the Gauss-Newton algorithm in Gallant 1987, 28-29
      for (iGN in 1:itmaxGN) {
        itersGN <- itersGN + 1;
        bprev2 <- bvec;
        
        ipmat <- probfunc(Y, Ypos, Xarray, tvec);
        presmat <- propmat - ipmat;
        loglik <- LogLik(Y,Ypos,ipmat,mvec)$LLu;
        if (print.level > 32 & iGN==1) print(paste("multinomMLE: -loglik initial:",loglik));
        score <-
          scorefunc(Ypos, nobs, tvunique, mvec, presmat, jacstack);
        hess2 <- hessianfunc(Ypos, nobs, tvunique, ipmat, mvec, jacstack) / nobs;
        posdef <- all(eigen(hess2, symmetric=TRUE, only.values=TRUE)$values > 0);
        gradient <- (score %*% rep(1,nobs)) / sqrt(sigma2GN) / nobs ;
        if (!posdef) {
          convflag <- FALSE;
          break;  ## quit if Hessian is not positive definite
        }
        bdiff <- c(solve(hess2, tol=.Machine$double.eps, LINPACK=TRUE) %*% gradient) ;  ## one Newton step
        ## print(bdiff);
        for (lambda in c(10:6)/10) {
          blambda <- bvec + lambda * bdiff;
          tlambda <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,blambda, ncats,tvars.total);
          plambda <- probfunc(Y, Ypos, Xarray, tlambda);
          logliklambda <- LogLik(Y,Ypos,plambda,mvec)$LLu;
          if (!is.na(logliklambda) && logliklambda < loglik) break;
        }
        if (is.na(logliklambda) || logliklambda > loglik) for (lambda in 2^(-(1:45))) {
          blambda <- bvec + lambda * bdiff;
          tlambda <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,blambda, ncats,tvars.total);
          plambda <- probfunc(Y, Ypos, Xarray, tlambda);
          logliklambda <- LogLik(Y,Ypos,plambda,mvec)$LLu;
          if (!is.na(logliklambda) && logliklambda < loglik) break;
        }
        if (logliklambda < loglik) bvec <- blambda;
        tvec <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,bvec, ncats,tvars.total);
        if (print.level > 32)  {
          cat("multinomMLE: ibvec: "); print(bvec)
        }
        convflag <- converged(bvec,bprev2) & converged2(logliklambda,loglik);
        ##     convflag <- convflag & all(abs(gradient) < 1e-9);
        if (convflag) break;
      }
      if (!posdef & print.level >= 0) {
        print("multinomMLE: Hessian is not positive definite");
      }
      if (print.level > 32 & posdef) {
        print(paste("multinomMLE: -loglik final: ",logliklambda));
      }
      LL2 <- ifelse(posdef, LogLik(Y,Ypos,plambda,mvec), NA);
      if (print.level > 32 & posdef) {
        print(paste("multinomMLE: -loglik:  unweighted,",LL2$LLu));
        print("multinomMLE: gradient:");  print(c(gradient));
        print("multinomMLE: bvec:");  print(bvec);
      }
      information <- hessianfunc(Ypos, nobs, tvunique, ipmat, mvec, jacstack);
      if (all(eigen(information, symmetric=TRUE, only.values=TRUE)$values > 0)) {
        formation <- solve(information, tol=.Machine$double.eps, LINPACK=TRUE);
      }
      else {
        formation <- NA;
      }
      return(
             list(coefficients=bvec, tvec=tvec, formation=formation, score=score,
                  LLvals=LL2, convflag=convflag, iters=itersGN, posdef=posdef) );
    }
  error <- 0;
  iters <- 0;
  for (i in 1:itmax) {
    iters <- iters + 1;
    bprev <- bvec;

    tvec <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,bvec, ncats,tvars.total);
    sigma2  <- sum(resfunc(Y,Ypos,Xarray,tvec)^2)/(nobs-tvunique);

    ## grouped multinomial:  estimate using Newton algorithm
    GNlist <- GNstep(bvec,sigma2);
    error <- ifelse(GNlist$posdef,0,32);  ## error == 32 if hessian not posdefinite
    if (print.level > 32)
      print(paste("multinomMLE: number of Newton iterations", GNlist$iters));

    bvec <- GNlist$coeff;
    if (converged(bvec,bprev)) break;
  }

  opg <- GNlist$score %*% t(GNlist$score) ;
  obsformation <- GNlist$formation ;
  rcovmat <- obsformation %*% opg %*% obsformation;

  if (print.level > 2) {  
    if (length(obsformation)==1 && obsformation==NA) {
      print(paste("multinomMLE: hessian determinant:",NA));
    }
    else {
      print(paste("multinomMLE: hessian determinant:",
        det(solve(obsformation, tol=.Machine$double.eps, LINPACK=TRUE))));
    }
    print(paste("multinomMLE: OPG determinant:", det(opg)));
  }

  ## table of returned error values (indicated values add to give total error)
  ## 0    no errors
  ## 32  Hessian not positive definite in the final Newton step

  se.opg.vec  <- sqrt(diag(ginv(opg/sigma2)));
  se.hes.vec <- sqrt(diag(obsformation));
  se.vec     <- sqrt(diag(rcovmat));

  se.opg <- xvec;
  se.hes <- xvec;
  se     <- xvec;    

  se.opg <- as.data.frame(mnl.xvec.mapping(forward=FALSE,xvec,se.opg,se.opg.vec,
                                           ncats,tvars.total));
  se.hes <- as.data.frame(mnl.xvec.mapping(forward=FALSE,xvec,se.hes,se.hes.vec,
                                           ncats,tvars.total));
  se <- as.data.frame(mnl.xvec.mapping(forward=FALSE,xvec,se,se.vec,
                                       ncats,tvars.total));
  tvec  <- as.data.frame(tvec)

  row.names(tvec) <- xvar.labels;
  names(tvec)     <- choice.labels;   

  row.names(se) <- xvar.labels;
  names(se)     <- choice.labels;

  row.names(se.opg) <- xvar.labels;
  names(se.opg)     <- choice.labels;

  row.names(se.hes) <- xvar.labels;
  names(se.hes)     <- choice.labels;

  return(
         list(coefficients=GNlist$tvec, coeffvec=GNlist$coeff, dispersion=sigma2,
              se=se, se.opg=se.opg, se.hes=se.hes,
              se.vec=se.vec, se.opg.vec=se.opg.vec,se.hes.vec=se.hes.vec,
              A=opg/sigma2, B=obsformation, covmat=rcovmat,
              iters=iters, error=error,
              GNlist=GNlist, sigma2=sigma2,
              Y=Y, Ypos=Ypos, fitted.prob=probfunc(Y, Ypos, Xarray, GNlist$tvec),
              jacstack=jacstack)
         )
}

#
#  multinomRob
#
#  Walter R. Mebane, Jr.
#  Cornell University
#  http://macht.arts.cornell.edu/wrm1/
#  wrm1@macht.arts.cornell.edu
#
#  Jasjeet Singh Sekhon 
#  Harvard University
#  http://jsekhon.fas.harvard.edu/
#  jsekhon@fas.harvard.edu
#
#  $Id: multinomRob.R,v 1.12 2004/02/19 02:13:11 wrm1 Exp $


######################################################
## multinomRob: genoud and Gauss-Newton estimation  ##
######################################################

multinomRob <-
  function(model,
           data,
           starting.values=NULL,
           equality=NULL,  # list of lists of parameter equality constraints
           genoud.parms   = NULL, ## can be NULL, single, or list
           print.level    = 0,
           iter = FALSE,   # should we iterate
           maxiter = 10,  # maximum number of iterations before we stop
           multinom.t=1,  # 0=no, 1=yes, 2=force should we do multinom-t for starting values
           multinom.t.df=NA #if set, the multivariate-t function is FORCED to use this DF 
           ) 
  {

    #check input
    if (print.level < 0)
      print.level  <- 0;

    if (iter!=TRUE & iter!=FALSE)
      {
        stop("multinomRob(): illegal input for variable iter")
      }    
    if (maxiter < 0)
      {
        stop("multinomRob(): illegal input for variable maxiter")
      }
    if (multinom.t!=0 & multinom.t!=1 & multinom.t!=2)
      {
        stop("multinomRob(): illegal input for variable multinom.t")
      }
    if (!is.na(multinom.t.df) & multinom.t.df <= 0)
      {
        stop("multinomRob(): illegal input for variable multinom.t.df")
      }    

    converged.test <- function(old,new,tol=1e-8) {
      abs((old-new)/old) < tol;
    }
    genoud.fun <- function(z) {
      fit.multinomial.C.lqd2(z,X,Y,Ypos,xvec,tvec,ncats,nvars,nvars.unique,obs,TotalY)
    }

    ## ----------------- BEGIN: PARSING / MODEL BUILDING AUTOMATION ------------------

    ## ---- DATA X, Y and Ypos ----
    data.all <- get.xy(model, data, print.level=print.level);
    Y <- data.all$Y;
    ncats  <- dim(Y)[2];
    obs    <- dim(Y)[1];
    Ypos <- data.all$ypos;  # element is TRUE if y>=0, FALSE if y<0
    Y[!Ypos] <- 0;  # set count values to be ignored to zero, for convenience later
    X <- data.all$X;
    nvars <- dim(X)[2];

    ## ---- LABELS AND CONTROLS/FLAGS ----
    xvec <- matrix(0,nrow=nvars,ncol=ncats)
    colnames(xvec) <-  choice.labels  <- data.all$ynames;

    xvar.labels <- rep( "", nvars)
    for (i in 1:ncats) {
      ## only look at named columns
      for (j in 1:nvars) {
        ## and keep only non-padding/zeroed names
        ssplit <- ifelse( i==1, "", "/")
        if (data.all$xlengths[i] >= j) {
          xvar.labels[j] <-
            paste( xvar.labels[j], ssplit, data.all$xnames[[i]][j], sep="")
          xvec[j,i] <- 1
        } else {
          xvar.labels[j] <- paste( xvar.labels[j], ssplit, "NA", sep="")
        }
      }
    }
    rownames(xvec) <- xvar.labels

    # parameter equality constraints
    if (!is.null(equality)) {
      eqlen <- length(equality);
      ieq <- 1;
      for (i in 1:eqlen) {
        ieq <- ieq + 1;
        xynames <- get.xynames(equality[[i]], data);
        nynames <- length(xynames$ynames);
        yidx <- match(xynames$ynames, data.all$ynames) ;
        if (any(is.na(yidx))) {
          if (print.level >= 0) {
            print("response used in equality constraints is not in the model");
            cat("responses in model:", data.all$ynames, "\n");
            cat("responses in equality:", xynames$ynames, "\n");
            print("terminating multinomRob")
          }
          return(NULL);
        }
        xidx <- list();
        xidxlen <- 0;
        for (j in 1:nynames) {
          xidx[[j]] <- match(xynames$xnames[[j]], data.all$xnames[[ yidx[j] ]]);
          if (any(is.na(xidx[[j]]))) {
            if (print.level >= 0) {
              print("regressor used in equality constraints is not in the model");
              cat("response:", xynames$ynames[j]);
              cat("regressors in model:", data.all$xnames[[ yidx[j] ]], "\n");
              cat("regressors in equality:", xynames$xnames, "\n");
              print("terminating multinomRob")
            }
            return(NULL);
          }
          xidxlen <- xidxlen + length(xidx[[j]]);
        }
        xyvmatrix <- matrix(0, 3, xidxlen);
        kk <- 0;
        for (j in 1:nynames) {
          for (k in 1:length(xidx[[j]])) {
            kk <- kk + 1;
            xyvmatrix[1,kk] <- xidx[[j]][k] ;
            xyvmatrix[2,kk] <- yidx[j] ;
            xyvmatrix[3,kk] <- xvec[xidx[[j]][k], yidx[j]] ;
          }
        }
        if (any(xyvmatrix[3,] != 1)) {  # need to consolidate equality constraints
          if (any(xyvmatrix[3,1] != xyvmatrix[3,])) {
            # not a subset of previous constraints
            xyvu <- unique(xyvmatrix[3,]);
            xyvu <- xyvu[xyvu>1];
            xyvmin <- min(xyvu);  # code to consolidate on
            for (k in 1:length(xyvu)) {
              xvec[xvec == xyvu[k]] <- xyvmin;
              xyvmatrix <- xyvmatrix[,xyvmatrix[3,] != xyvu[k]];
            }
            if (dim(xyvmatrix)[2] > 0) { # parameters not previously constrained
              for (k in 1:(dim(xyvmatrix)[2])) {
                xvec[xyvmatrix[1,k], xyvmatrix[2,k]] <- xyvmin;
              }
            }
          }
        }
        else { # consolidation not needed
          for (k in 1:(dim(xyvmatrix)[2])) {
            xvec[xyvmatrix[1,k], xyvmatrix[2,k]] <- ieq;
          }
        }
      }      
      if (print.level > 0) {
        cat("\nEquality constraints among parameters (after consolidation):\n");
        nidxvals <- length(idxvals <- sort(unique(xvec[xvec>1])));
        for (i in 1:nidxvals) {
          cat("Equality constrained set", i, "\n");
          for (j in 1:nvars)  for (k in 1:ncats) {
            if (xvec[j,k]==idxvals[i]) {
              cat("outcome", data.all$ynames[k], "regressor",
                  data.all$xnames[[k]][j], "\n");
            }
          }
        }
      }
    }
    
    if (print.level > 0) {
      cat("\nYour Model (xvec):\n")
      print(xvec)
    }

    TotalY <- apply( Y, 1, sum)
    tvec  <- xvec;
    nvars.unique <- sum(xvec == 1) + length(unique(xvec[xvec>1]));
    nvars.total  <- nvars

    ## ------------------- END: PARSING / MODEL BUILDING AUTOMATION ------------------

    #let's create jacstack
    jacstack  <- jacstack.function(X,nvars.unique,xvec)

    # check for regressors with a distinct value at only one observation
    jsingle <- jacstack.singles(jacstack);
    if (any(jsingle) & print.level > 0) {
      cat("\n");
      print("multinomRob:  WARNING.  Limited regressor variation...")
      print("WARNING.  ... A regressor has a distinct value for only one observation.")
      print("WARNING.  ... I'm using a modified estimation algorithm (i.e., preventing LQD")
      print("WARNING.  ... from modifying starting values for the affected parameters).")

      xsingle <- tvec != tvec;
      xsingle <- mnl.xvec.mapping(forward=FALSE,xvec,xsingle,jsingle,ncats,nvars);
      print("WARNING.  ... Affected parameters are TRUE in the following table.")
      cat("\n");
      print(xsingle);
      cat("\n");
    }

    #load up genoudParms
    if(is.null(genoud.parms))
      genoud.parms  <- list(Domains=NULL);
    genoud.parms  <- genoudParms(genoud.parms)

    # initialize with really bad fit values
    mnl.fit  <- 9999999999;
    multinomT.fit  <- 9999999999;
    starting.fit   <- 9999999999;
    multinomT.foo <- NULL
    if (is.null(starting.values)) {
      starting.values <- vector(mode="numeric", length=nvars.unique);      
      if (print.level > 0)
        cat("\n\nmultinomRob(): Grouped MNL Estimation\n");
      Yp  <- Y/TotalY
      Yp[,ncats]  <- 1-apply(as.matrix(Yp[,1:(ncats-1)]),1,sum);
      mnl1 <- multinomMLE(Y=Y, Ypos=Ypos, Xarray=X, xvec=xvec, jacstack=jacstack,
                          xvar.labels=xvar.labels, choice.labels=choice.labels,
                          print.level=print.level);

      starting.values <- mnl1$coeffvec
      mnl.fit <-
        fit.multinomial.C.lqd2(starting.values,X,Y,Ypos,xvec,tvec,ncats,nvars,
                               nvars.unique,obs,TotalY);
      
      if (print.level > 0) {
        cat("MNL LQD Fit:",mnl.fit,"\n");        
        cat("MNL Estimates:\n");
        print(mnl1$coefficients);
        cat("\n");
        cat("MNL SEs:\n");
        print(mnl1$se);        
        cat("\n");
      }

      starting.values.user  <- starting.values
      tvec    <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,starting.values,
                                  ncats,nvars)
      starting.fit <-
         fit.multinomial.C.lqd2(starting.values,X,Y,Ypos,xvec,tvec,ncats,nvars,
                                nvars.unique,obs,TotalY)

      use.mnl  <- TRUE;
      if (multinom.t > 0 & all(Ypos)) {
        if (print.level > 0) {        
          cat("\n\nmultinomRob(): Calculating multinomial-t starting values.\n");
        }

        Sres.raw <- res.std(Y,TotalY,mnl1$fitted.prob);
        multinomT.startvalues  <- list()
        multinomT.startvalues$beta  <- starting.values;
        multinomT.startvalues$Omega <- var(Sres.raw);
        #multinomT.startvalues$Omega <- var(mnl1$res[,1:(ncats-1)])/(obs-nvars.unique);
        multinomT.startvalues$df    <- 10;
        multinomT.foo  <- multinomT(Yp=Yp, Xarray=X, xvec=xvec, jacstack=jacstack,
                                    start=multinomT.startvalues, nobsvec=TotalY,
                                    fixed.df=multinom.t.df);

        for (ii in 0:6) {
          if(is.null(multinomT.foo$se$beta[1]) & is.na(multinom.t.df)) {
            if (print.level > 2)
              cat("Fixing Bad Estimate:",ii,"\n")
            multinomT.foo  <- multinomT(Yp=Yp, Xarray=X, xvec=xvec, jacstack=jacstack,
                                        start=multinomT.startvalues, nobsvec=TotalY,
                                        fixed.df=(10^ii-.5));
            if (print.level > 2)
              cat("DF: ",multinomT.foo$par$df,"\n");    
            
            multinomT.startvalues$beta  <- multinomT.foo$par$beta
            multinomT.startvalues$Omega <- multinomT.foo$par$Omega
            multinomT.startvalues$df    <- multinomT.foo$par$df
            multinomT.foo  <- multinomT(Yp=Yp, Xarray=X, xvec=xvec, jacstack=jacstack,
                                       start=multinomT.startvalues, nobsvec=TotalY);
          } else if (is.null(multinomT.foo$se$beta[1]) & !is.na(multinom.t.df)) {
            if (print.level > 1)
              {
                cat("WARNING: Multinom-T SEs are null, but we are using multinomT point estimates for starting values anyways because you explicitly sepecifed a DF\n")
              }
            use.mnl  <- FALSE;
            starting.values <- multinomT.foo$par$beta;
            multinomT.fit <-
              fit.multinomial.C.lqd2(starting.values,X,Y,Ypos,xvec,tvec,ncats,nvars,
                                     nvars.unique,obs,TotalY);
            if (print.level > 1)
              {
                cat("Multinom-T LQD Fit (step ",ii,"):",multinomT.fit,"\n")  
                beta.print <-
                  mnl.xvec.mapping(forward=FALSE,xvec,tvec,multinomT.foo$par$beta,ncats,nvars)
                cat("Multinom-T Beta Estimates (step ",ii,"):\n")
                print( beta.print )
                cat("Multinom-T SEs are not defined\n")
              }
            break;
          } else  {
            use.mnl  <- FALSE;
            starting.values <- multinomT.foo$par$beta;
            multinomT.fit <-
              fit.multinomial.C.lqd2(starting.values,X,Y,Ypos,xvec,tvec,ncats,nvars,
                                     nvars.unique,obs,TotalY);
            if (print.level > 1)
              {
                cat("Multinom-T LQD Fit (step ",ii,"):",multinomT.fit,"\n")  
                beta.print <-
                  mnl.xvec.mapping(forward=FALSE,xvec,tvec,multinomT.foo$par$beta,ncats,nvars)
                se.print <-
                  mnl.xvec.mapping(forward=FALSE,xvec,tvec,multinomT.foo$se$beta,ncats,nvars)
                cat("Multinom-T Beta Estimates (step ",ii,"):\n")
                print( beta.print )
                cat("Multinom-T Beta SEs (step ",ii,"):\n")
                print( se.print )
                cat("Multinom-T Omega Estimates (step ",ii,"):\n")
                print( multinomT.foo$par$Omega )
                cat("Multinom-T DF (step ",ii,"):", multinomT.foo$par$df,"\n\n")
              }
            break;
          }
        }#ii
      }#end of multinom.t

      if (use.mnl==FALSE & !is.null(starting.values) & multinom.t < 2) {
        if (multinomT.fit >= mnl.fit)
          use.mnl <- TRUE;
      }

      if(use.mnl==TRUE) {
        if (print.level > 0)
          cat("multinomRob(): Using grouped MNL estimates as starting values.\n")
        
        starting.values <- vector(mode="numeric", length=nvars.unique);
        starting.values <-
          mnl.xvec.mapping(forward=TRUE,xvec,mnl1$coefficients,starting.values,
                           ncats,nvars);
      } else {
        if (print.level > 0) {
          cat("multinomRob(): Using multinomial-t estimates as starting values.\n")
          cat("Multinom-T LQD Fit:",multinomT.fit,"\n")  
          beta.print <-
            mnl.xvec.mapping(forward=FALSE,xvec,tvec,multinomT.foo$par$beta,ncats,nvars)
          cat("Multinom-T Estimates:\n")
          print( beta.print )
          cat("Multinom-T DF:", multinomT.foo$par$df,"\n\n")
        }
        starting.values <- multinomT.foo$par$beta        
      }
    }#is.null(starting.values)
      
    tvec <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,starting.values,ncats,nvars);

    fit.starting <-
      fit.multinomial.C.lqd2(starting.values,X,Y,Ypos,xvec,tvec,ncats,nvars,
                             nvars.unique,obs,TotalY)

    if (print.level > 2) {
      cat("multinomRob(): Starting Values \n")
      print(tvec)
      cat("multinomRob(): starting fit =",fit.starting,"\n");
    }
        
    genoud.out  <- genoudRob(genoud.fun,
                          nvars=nvars.unique,
                          starting.values = starting.values,
                          as.list(genoud.parms))
      
    s0    <- genoud.out$value;
    lqd.beta.vector  <- genoud.out$par;
    tvec  <- xvec;
    tvec  <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,lqd.beta.vector,ncats,nvars);

    if (print.level>0) {    
      cat("\n");
      cat("\nLQD Results:\n");
      print(tvec);
      cat("\nLQD sigma:",s0,"\n\n")
    }
    
    #############################################
    ## TANH                                     #
    #############################################
    
    cat("\n(multinomTanh):\n");
      
    # modify LQD results for single unique value regressors
    tvec  <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,
                              ifelse(jsingle, starting.values, lqd.beta.vector),
                              ncats,nvars);

    mout <- multinomTanh(Y=Y, Ypos=Ypos, X=X, jacstack=jacstack,
                         xvec=xvec, tvec=as.data.frame(tvec), 
                         pop=TotalY,s2=s0^2,
                         xvar.labels=xvar.labels, choice.labels,
                         print.level=print.level);
    error <- mout$mtanh$error;        
    if (mout$mtanh$error==0) {
      if (print.level > 0) {
        cat("Tanh Estimates\n")
        print(mout$mtanh$coef)
        cat("\nTanh Sandwich SEs\n");
        print(mout$mtanh$se);
        cat("\n")
        cat("TANH sigma:",sqrt(mout$mtanh$tanhsigma2),"\n\n")        
      } #end print      
    } else {
      if (print.level >= 0) {
        cat("WARNING: Tanh returned an error code:",mout$mtanh$error,"\n");
        cat("WARNING: Optimization not complete.  Starting values are probably not good enough. \n")
        if(iter) {
          cat("WARNING: Hopefully we will fix this problem during subsequent iterations.\n")
        } else {
          cat("WARNING: Rerun with new different starting value or with the iteration option left on\n")
        }
      }
    }

    # If the fit value gets worse, we stop.  Which is *NOT*
    # what the reference file does.  The reference file
    # keeps going (even if the new value is worse until there is
    # stability.
    s.count  <- 0;
    if (iter)
      {
        s0.old  <- s0;
        s0.tanh.old  <- sqrt(mout$mtanh$tanhsigma2);
        tolerance  <- genoud.parms$solution.tolerance;

        #let's shrink our domains by a half for the iteration stuff because we are
        # going to assume that the tanh starting values are not that bad.
        genoud.parms$scale.domains  <- genoud.parms$scale.domains/2;

        mout.old  <- mout;
      }#iter
    while(iter)
      {
        s.count  <- s.count+1;
        if (print.level > 0) {
          cat("\n**********************************\n\n");
          cat("Iteration:",s.count,"\n")
        }

        if (!any(is.na(mout.old$mtanh$coeffvec))) {
          LQDfit.tanhcoefs  <-
            fit.multinomial.C.lqd2(mout.old$mtanh$coeffvec,X,Y,Ypos,xvec,tvec,
                                   ncats,nvars,nvars.unique,obs,TotalY);
        } else {
          LQDfit.tanhcoefs  <- 9999999999;
        }

        best  <- order(c(starting.fit, mnl.fit, multinomT.fit, LQDfit.tanhcoefs))[1]
        if (s.count > 1 && !any(is.na(mout.old$mtanh$coeffvec))) {
          if (print.level > 0)
            cat("using LQD sigma provided by Tanh.\n")
          starting.values  <- mout.old$mtanh$coeffvec
        } else if (best==1) {
          if (print.level > 0)
            cat("using best non-genoud LQD sigma. Provided by user starting values.\n")
          starting.values  <- starting.values.user
        } else if (best==2) {
          if (print.level > 0)
            cat("using best non-genoud LQD sigma. Provided by ML multinomial.\n")
          starting.values  <- mnl1$coeffvec
        }  else if (best==3) {
          if (print.level > 0)
            cat("using best non-genoud LQD sigma. Provided by multinomial-t.\n")
          starting.values  <- multinomT.foo$par$beta;
        } else if (best==4) {
          if (print.level > 0)
            cat("using best non-genoud LQD sigma. Provided by Tanh.\n")
          starting.values  <- mout.old$mtanh$coeffvec
        }

        genoud.out  <- genoudRob(genoud.fun,
                                 nvars=nvars.unique,
                                 starting.values = starting.values,
                                 as.list(genoud.parms))
        
        s0   <- genoud.out$value;
        lqd.beta.vector  <- genoud.out$par;
        tvec  <- xvec;
        tvec  <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,lqd.beta.vector,
                                  ncats,nvars);
        
        if (print.level>0) {    
          cat("\n(iter",s.count,"): LQD Results:\n");
          print(tvec);
          cat("\n(iter",s.count,"): LQD sigma:",s0,"\n\n")
        }

        if (s0 > (s0.old+tolerance))
            {
              if (print.level > 0)
                {
                  cat("WARNING: New fit is worse.  Stopping.\n");
                  cat("It may be a good idea to restart with a larger GENOUD population.\n");
                } #end of print.level
              error  <- -1;
              iter  <- FALSE;              
              break;
            }

        if (!is.na(s0.tanh.old) && converged.test(s0.old, s0, tolerance))
          {
              if (print.level > 0)
                {
                  cat("\n**********************************\n")
                  cat("CONVERGED\n");

                  if (print.level > 0) {
                    cat("(Converged): Tanh Estimates\n")
                    print(mout$mtanh$coef)
                    cat("\n(Converged): Tanh Sandwich SEs\n");
                    print(mout$mtanh$se);
                    cat("\n")
                    cat("(Converged): TANH sigma:",sqrt(mout$mtanh$tanhsigma2),"\n\n")
                    s0.tanh  <- sqrt(mout$mtanh$tanhsigma2);
                  } #end print                                    
                } #end of print.level
              error  <- 0;
              iter  <- FALSE;              
              break;            
          } else {
            #new s0 is smaller
        
            # modify LQD results for single unique value regressors
            tvec  <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,
                                 ifelse(jsingle, starting.values, lqd.beta.vector),
                                 ncats,nvars);

            mout <- multinomTanh(Y=Y, Ypos=Ypos, X=X, jacstack=jacstack,
                                 xvec=xvec, tvec=as.data.frame(tvec), 
                                 pop=TotalY,s2=s0^2,
                                 xvar.labels=xvar.labels, choice.labels,
                                 print.level=print.level);
            
            error  <- mout$mtanh$error;        
            if (mout$mtanh$error==0)
              {
                if (print.level > 0) {
                  cat("(iter",s.count,"): Tanh Estimates\n")
                  print(mout$mtanh$coef)
                  cat("\n(iter",s.count,"): Tanh Sandwich SEs\n");
                  print(mout$mtanh$se);
                  cat("\n")
                  cat("(iter",s.count,"): TANH sigma:",sqrt(mout$mtanh$tanhsigma2),"\n\n")
                  s0.tanh  <- sqrt(mout$mtanh$tanhsigma2);

                  mout.old  <- mout;                  
                  s0.tanh.old  <- s0.tanh;
                } #end print                  
              } else {
                if (print.level > 0)
                  cat("(iter",s.count,"): tanh returns an error:",mout$mtanh$error,"\n");
              }
          } #end of else
            
        if (s.count==maxiter)
          {
            if (print.level > 0) 
              cat("\nWARNING: Maximum number of iterations reached.  Results have NOT converged\n\n");

            #iter failed
            error  <- -1;
            iter  <- FALSE;
            break;
          }
        if (s0 < s0.old)
          s0.old  <- s0;        
      }#iter check

    if (!exists("mnl1"))
      mnl1  <- NULL;
    if (!exists("multinomT.foo"))
      multinomT.foo  <- NULL;

    #Rotated Residuals.  This will result in a relatively easy to interpret
    #vector of residuals.  But this residual vector will NOT be a
    #consistent set of ortho residuals.
    residuals.rotate  <- matrix(nrow=obs,ncol=ncats)
    for (ii in 1:ncats)
      {
        tindx  <- 1:ncats
        tindx[1]  <- ii;
        tindx[ii] <- 1;
        
        YTmp     <- Y[,c(tindx)]
        YposTmp  <- Ypos[,c(tindx)]
        XTmp     <- X[,,c(tindx)]
        jacstackTmp  <- jacstack[,,c(tindx)]
        tvec  <- as.matrix(mout$mtanh$coefficients[,c(tindx)])
        
        foo  <- permute(Y=YTmp, Ypos=YposTmp, X=XTmp,
                        jacstack=jacstackTmp, tvec=tvec,
                        pop=TotalY, sigma=sqrt(mout$mtanh$disp),
                        weight=mout$mtanh$w)
        
        residuals.rotate[,ii]  <- foo$student[,1];
      }  #end of ii loop
    
    z  <- list(coefficients=mout$mtanh$coefficients,
               se=mout$mtanh$se,
               LQDsigma2=mout$mtanh$dispersion,
               TANHsigma2=mout$mtanh$tanhsigma2,
               weights=mout$weights,
               Hdiag=mout$Hdiag,
               prob=mout$mtanh$prob,
               residuals.rotate=residuals.rotate,
               residuals.student=mout$cr$student,
               residuals.standard=mout$cr$standard,
               mnl=mnl1,
               multinomT  = multinomT.foo,
               genoud = genoud.out,
               mtanh = mout$mtanh,
               error = error,
               iter = s.count); #the iter number at the solution
    class(z)  <- "multinomRob"
    return(z)
  }


#mimicking the definition of summary.default
#summary.default <-
#    function(object, ..., digits = max(3, getOption("digits") - 3))
summary.multinomRob <- function(object, ..., digits=3, weights=FALSE)
{
  
  
  if (class(object) != "multinomRob") {
    warning("Object not of class 'multinomRob'")
    return(NULL)
  } 

  out.mtanh <- list()

  if (!is.null(object$mtanh) && is.list(object$mtanh) && object$mtanh$error == 0) {
    ncats <- NCOL(object$mtanh$coefficients);
    nx <- NROW(object$mtanh$coefficients);
    # find length of longest regressor label
    labmax <- 0;
    for (i in 1:ncats) {
      for (j in 1:nx) {
        labn <- nchar(strsplit(row.names(object$mtanh$se)[j], "/")[[1]][i]) ;
        if (labmax < labn) labmax <- labn;
      }
    }
    # generate a blank variable with as many spaces as the longest label
    spc <- " ";
    blank <- "";
    for (i in 1:labmax) blank <- paste(spc, blank, sep="");

    for (i in 1:ncats) {
      tmp <- as.data.frame(
               list("Est" = object$mtanh$coefficients[,i],
                    "SE Sand" = object$mtanh$se[,i],
#                    "SE OPG"  = object$mtanh$se.opg[,i],
#                    "SE Hess" = object$mtanh$se.hes[,i],
                    "t-val Sand" = object$mtanh$coefficients[,i]/object$mtanh$se[,i]))
      tmp <- as.matrix(tmp);
      rn <- rep(blank, nx);
      for (j in 1:nx) {
        rnj <- strsplit(row.names(object$mtanh$se)[j], "/")[[1]][i] ;
        substr(rn[j], 1, nchar(rnj)) <- rnj;
      }
      dimnames(tmp)[[1]] <- rn;
      out.mtanh[[i]] <- tmp
      choice.lables  <- labels(object$coefficients)[[2]]
      cat("\nChoice",i,":",choice.lables[i],"Estimates and SE:\n")
      print(signif(tmp,digits=digits))
      cat("\n")
    }

    cat("\n");
    cat("LQD sigma:",sqrt(object$LQDsigma2),"\n")
    cat("TANH sigma:",sqrt(object$TANHsigma2),"\n\n")
        
    indx  <- object$weights[,2:ncats]==0;
    indx[is.na(indx)] <- FALSE;
    if (ncats==2) indx <- matrix(indx, length(indx), 1);
    w0.obs    <- sum(apply(indx,1,any));
    cat("Number of Observations:",NROW(object$prob),"\n")
    cat("Number of observations with at least one zero weight:",w0.obs,"\n")
    w0  <- sum(indx)
    cat("Number of zero weights:",w0,"\n")

    if (weights) {
      cat("TANH: weights\n");
      print(object$weights);
    }
  }
  else {
    print("error encountered in tanh result object");
  }
} #end of summary.multinomRob()


#mimicking the definition of plot 
plot.multinomRob  <- function(x, ...)
  {
    #this plots the residuals
    if (class(x) != "multinomRob") {
      warning("Object not of class 'multinomRob'")
      return(NULL)
    } 

    out.mtanh <- list()

    if (!is.null(x$mtanh) && is.list(x$mtanh) && x$mtanh$error == 0)
      {
        ncats <- NCOL(x$mtanh$coefficients);        
        choice.lables  <- labels(x$coefficients)[[2]]

        ask.start  <- par("ask",no.readonly=TRUE)

        for (i in 1:(ncats-1))
          {
            plot(x$residuals.student[,i],xlab="Observation", ylab="Residual", ...)
            title(main=paste(choice.lables[i],"\nStudentized Residuals"));

            if (i==1)
              par(ask=TRUE);
          } #end for
        par(ask=ask.start)
      }
    else {
      print("error encountered in tanh result object");
    }
    #end if
  }#end of plot.multinomRob


#Function to create ortho-residuals (with base correction) for the
#other choices.  This will result in a relatively easy to interpret
#vector of residuals.  But this residual vector will NOT be a
#consistent set of ortho residuals.
#
#Modeled on
#lapo:~/xchg/election/R/multinomial/FL2/base4.permute4.origtanh2.R
#permute.newtanh4.R
#
permute  <- function(Y, Ypos, Xarray, jacstack, tvec, pop, sigma, weight)
  {
#tvec: tanh estimated values    
    #from multinomTanh
    nobs  <- dim(Y)[1]
    ncats <- dim(Y)[2]
    Hdiag <- robustified.leverage(tvec, Y, Ypos, Xarray, pop, ifelse(weight >0,1,0),jacstack);
    w.Hdiag <- as.data.frame(matrix(c(as.vector(1:nobs),
                                     signif(weight), signif(Hdiag)),ncol=(ncats-1)+(ncats-1)+1));
#    names(w.Hdiag) <-
#      c("name",paste("weights:",choice.labels[1:ncats-1],sep=""),
#        paste("Hdiag:",choice.labels[1:ncats-1],sep=""));
    #cat("mtanh: weights, Hdiag (by choices)\n");

    cr <- fn.region.results(tvec, Y, Ypos, Xarray, pop, sigma, Hdiag);    
    return(list(pred=cr$pred, student=cr$student, standard=cr$standard, Hdiag=Hdiag))
  } #end of permute
#
#  multinomRob
#
#  Walter R. Mebane, Jr.
#  Cornell University
#  http://macht.arts.cornell.edu/wrm1/
#  wrm1@macht.arts.cornell.edu
#
#  Jasjeet Singh Sekhon 
#  Harvard University
#  http://jsekhon.fas.harvard.edu/
#  jsekhon@fas.harvard.edu
#
#  $Id: multinomT.R,v 1.6 2004/02/19 02:13:11 wrm1 Exp $
#
#

##   Yp:  matrix of (overdispersed and contaminated) multinomial proportions
##   Xarray:  array of regressors,
##      dim(Xarray) = c(n observations, n parameters, n categories)
##   xvec:  vector to indicate all the coefficient parameters in the model
##      (parms by ncats):
##      It has a 1 for an estimated parameter and a 0 otherwize.
##      example:
##      > xvec
##           [,1] [,2] [,3] [,4] [,5]
##      [1,]    1    1    1    1    0
##      [2,]    1    1    1    1    0
##      [3,]    1    1    1    1    0
##      [4,]    1    1    1    1    0
##   tvec: parms by ncats matrix (matrix with LQD estimates):
##      example:
##      > tvec
##                          Buchanan        Nader     Gore     Bush Other
##      int               -0.1641034    1.0735560 3.363641 4.151853     0
##      p(r,dg,d,r)96      2.4413780    0.3269827 3.207104 1.676676     0
##      p(cr,g,cd,cr)RV00 15.5333800 1149.4130000 2.039766 1.761392     0
##      pCuban            -6.4083750    0.3546630 1.966287 2.795598     0
##   jacstack:  array of regressors,
##      dim(Xarray) = c(n observations, n UNIQUE parameters, n categories)

multinomT  <- function(Yp, Xarray, xvec, jacstack,
                       start=NA, nobsvec, fixed.df = NA)
  {
    #Y (raw-Y) is assumed to be proportions,
    #the last category to be the contrast
    
    obs    <- dim(Yp)[1];
    cats   <- dim(Yp)[2];
    mcats  <- cats-1;
    nvars  <- dim(Xarray)[2];

    if (any(xvec[,cats] != 0)) {
      stop("(multinomT): invalid specification of Xarray (regressors not allowed for last category");
    }
    smdim <- dim(jacstack);
    stack.index <- matrix(FALSE, smdim[2], smdim[3]);
    for (i in 1:smdim[2]) for (j in 1:smdim[3]) {
      stack.index[i,j] <- !all(jacstack[,i,j] == 0);
    }
    if (sum(stack.index) != sum(xvec != 0)) {
      print("multinomT:  xvec is:"); print(xvec);
      print("multinomT:  stack.index is:"); print(stack.index);
      stop("(multinomT):  jacstack structure check failed");
    }
    kY  <- matrix(nrow=obs,ncol=mcats)
    indx1  <- Yp==0;
    indx1vec  <- apply(indx1,1,sum) > 0;
    if (sum(indx1) > 0)
      {
        cat("multinomT:  Need to remove 0 in multinomT transformation\n");
      }
    Yp[indx1]  <- .5/nobsvec[indx1vec];
    kY  <- log(Yp[,1:mcats]/Yp[,cats])
#    indx  <- is.infinite(kY);

    if (is.na(start[1]))
      {
        mt.obj  <- mt.mle(Xarray=Xarray, xvec=xvec, jacstack=jacstack, y=kY,
                          stack.index=stack.index, fixed.df=fixed.df);
      } else {
        mt.obj  <- mt.mle(Xarray=Xarray, xvec=xvec, jacstack=jacstack, y=kY,
                          stack.index=stack.index, start=start, fixed.df=fixed.df);        
      }

    tvec <- mnl.xvec.mapping(forward=FALSE,xvec,xvec, mt.obj$dp$beta, cats, nvars);
    pred <- mnl.probfunc(Yp, Yp==Yp, Xarray, tvec)
    
    return( list(call=mt.obj$call, logL=mt.obj$logL, deviance=mt.obj$deviance,
           par=mt.obj$dp, se=mt.obj$se, optim=mt.obj$optim, pred=pred))
  }#end of multinomT


#
#this version includes analytical gradients, no bounds on DF, other than df>0
#

mt.mle  <- function (Xarray, xvec, jacstack, y, stack.index,
          start = NA, freq =NA, fixed.df = NA, trace = FALSE,
          method = "BFGS", control = list(maxit = 600,trace=0, tol=1e-100)) 
{
  nvars  <- dim(Xarray)[2];
  Diag <- function(x) diag(x, nrow = length(x), ncol = length(x))
  y.name <- deparse(substitute(y))
  y.names <- dimnames(y)[[2]]
  y <- as.matrix(y)
  if (missing(freq) | is.na(freq)) {
#  if (missing(freq)) {    
    freq <- rep(1, nrow(y))
  }
#  x.names <- dimnames(mX)[[2]]

  d <- ncol(y)
  n <- sum(freq)
  m <- sum(xvec == 1) + length(unique(xvec[xvec>1]));

  if (is.na(start[1]))
    {
      cat("mt.mle:  I don't know how to generate starting values in the general case\n");
      stop();
    } #end if
  beta <- start$beta
  Omega <- start$Omega
  if (!is.na(fixed.df)) start$df <- fixed.df;
  df <- start$df

  Oinv <- solve(Omega, tol=1e-100)
  Oinv <- (Oinv + t(Oinv))/2
  upper <- chol(Oinv)
  D <- diag(upper)
  A <- upper/D
  D <- D^2
  if (d > 1)
    {
      param <- c(beta, -0.5 * log(D), A[!lower.tri(A, diag = TRUE)])
    } else {
      param <- c(beta, -0.5 * log(D))
    }
  if (is.na(fixed.df)) 
    param <- c(param, log(df))
  opt <- optim(param, fn = mt.dev, method = method, 
               control = control, hessian = TRUE,
               Xarray=Xarray, xvec=xvec, jacstack=jacstack, y = y,
               stack.index = stack.index, nvars = nvars, freq = freq,
               trace = trace, fixed.df = fixed.df)
  dev <- opt$value
  param <- opt$par
  if (trace) {
    cat("mt.mle:  Message from optimization routine:", opt$message, 
        "\n")
    cat("mt.mle:  deviance:", dev, "\n")
  }
  beta <- param[1:m] ;
  D <- exp(-2 * param[(m + 1):(m + d)])
  if (d > 1) {
    A <- diag(d)
    A[!lower.tri(A, diag = TRUE)] <-
      param[(m + d + 1):(m + d + d * (d - 1)/2)]
    i0 <- m + d + d * (d - 1)/2
  } else {
    i0 <- m + 1
    A <- as.matrix(1)
  }
  if (is.na(fixed.df))
    {
        df <- exp(param[i0 + 1])
      } else {
        df <- fixed.df
      }
  Ainv <- backsolve(A, diag(d))
  Omega <- Ainv %*% Diag(1/D) %*% t(Ainv)
  omega <- sqrt(diag(Omega))
#  dimnames(beta) <- list(x.names, y.names)
  dimnames(Omega) <- list(y.names, y.names)
  info <- opt$hessian/2
  if (all(is.finite(info))) {
    qr.info <- qr(info)
    info.ok <- (qr.info$rank == length(param))
  } else {
      info.ok <- FALSE
    }
  if (info.ok) {
    se2 <- diag(solve(qr.info))
    if (min(se2) < 0)
      {
        se <- NA
          } else {
            se <- sqrt(se2)
            se.beta <- se[1:m] ;
#            dimnames(se.beta)[2] <- list(y.names)
#            dimnames(se.beta)[1] <- list(x.names)
            se.df <- df * se[i0 + 1]
            se <- list(beta = se.beta, df = se.df, 
                       info = info)
          }
  } else {
    se <- NA
  }
  dp <- list(beta = beta, Omega = Omega, df = df)
  list(call = match.call(), logL = -0.5 * dev, deviance = dev, 
       dp = dp, se = se, optim = opt)
}


mt.dev  <- function (param, Xarray, xvec, jacstack, y, stack.index, nvars, freq,
    fixed.df = NA, trace = FALSE) 
{
    Diag <- function(x) diag(x, nrow = length(x), ncol = length(x))
    d <- ncol(y)
    n <- sum(freq)
    m <- sum(xvec == 1) + length(unique(xvec[xvec>1]));
    beta <- param[1:m];
    D <- exp(-2 * param[(m + 1):(m + d)])
    if (d > 1) {
        A <- diag(d)
        A[!lower.tri(A, diag = TRUE)] <-
          param[(m + d + 1):(m + d + d * (d - 1)/2)]
        i0 <- m + d + d * (d - 1)/2
    } else {
        i0 <- m + 1
        A <- as.matrix(1)
    }
    eta <- rep(0,d);
    if (is.na(fixed.df))
      {
        df <- exp(param[i0 + 1])
      } else {
        df <- fixed.df
      }
    Oinv <- t(A) %*% Diag(D) %*% A
#    cat("beta:\n");
#    print(as.matrix(beta));
    tvec <- mnl.xvec.mapping(forward=FALSE, xvec, xvec, beta, d+1, nvars);
    ylinpred <- y;
    for (j in 1:d) {
      ylinpred[,j] <- Xarray[,,j] %*% tvec[,j];
    }
    u <- y - ylinpred ;
#    cat("u:\n")
#    print(u)
#    cat("boo1\n");
    Q <- apply((u %*% Oinv) * u, 1, sum)
    L <- as.vector(u %*% eta)
    logDet <- sum(log(df * pi/D))
    dev <- (n * (2 * lgamma(df/2) + logDet - 2 * lgamma((df + 
        d)/2)) + (df + d) * sum(freq * log(1 + Q/df)) - 2 * sum(freq * 
        log(2 * pt(L * sqrt((df + d)/(Q + df)), df + d))))
    if (trace) 
        cat("mt.dev: ", dev, "\n")
    dev
}


mt.dev.grad  <- function (param, Xarray, xvec, jacstack, y, stack.index, nvars, freq,
                          fixed.df = NA, trace = FALSE) 
{
    Diag <- function(x) diag(x, nrow = length(x), ncol = length(x))
    d <- ncol(y)
    n <- sum(freq)
    m <- sum(xvec == 1) + length(unique(xvec[xvec>1]));
    nvarsunique <- dim(jacstack)[2];
    beta <- param[1:m];    
    D <- exp(-2 * param[(m + 1):(m + d)])
    if (d > 1) {
        A <- diag(d)
        A[!lower.tri(A, diag = TRUE)] <-
          param[(m + d + 1):(m + d + d * (d - 1)/2)]
        i0 <- m + d + d * (d - 1)/2
    }
    else {
        i0 <- m + d
        A <- as.matrix(1)
    }
    eta <- rep(0,d);
    if (is.na(fixed.df)) 
        df <- exp(param[i0 + 1])
    else df <- fixed.df
    tA <- t(A)
    Oinv <- tA %*% Diag(D) %*% A
    tvec <- mnl.xvec.mapping(forward=FALSE, xvec, xvec, beta, d+1, nvars);
    ylinpred <- y;
    for (j in 1:d) {
      ylinpred[,j] <- Xarray[,,j] %*% tvec[,j];
    }
    u <- y - ylinpred ;
    Q <- as.vector(apply((u %*% Oinv) * u, 1, sum))
    L <- as.vector(u %*% eta)
    t. <- L * sqrt((df + d)/(Q + df))
    dlogft <- -(df + d)/(2 * df * (1 + Q/df))
    dt.dL <- sqrt((df + d)/(Q + df))
    dt.dQ <- (-0.5) * L * sqrt(df + d)/(Q + df)^1.5
    T. <- pt(t., df + d)
    dlogT. <- dt(t., df + d)/T.
    u.freq <- u * freq
#    foo1 <- foo2 <- foo3 <- matrix(0, nvarsunique, d+1);
    fooA <- matrix(0, nvarsunique, d);
    for (j in 1:d) {
      fooA[,j] <- t(jacstack[,,j]) %*% (u.freq * (dlogft + dlogT. * dt.dQ));
#      foo1[,j] <- t(jacstack[,,j]) %*% (u.freq * dlogft);
#      foo3[,j] <- t(jacstack[,,j]) %*% (dlogT. * dt.dQ * u.freq);
#      foo2[,j] <- t(jacstack[,,j]) %*% (dlogT. * dt.dL * freq);
    }
#    foo1 <- t(mX) %*% (u.freq * dlogft);
#    foo2 <- t(mX) %*% (dlogT. * dt.dL * freq);
#    foo3 <- t(mX) %*% (dlogT. * dt.dQ * u.freq);
#    Dbeta <- (-2 * (foo1 + foo3) %*% Oinv) - outer(as.vector(foo2), eta)
    Dbeta <- -2 * fooA %*% Oinv
    if (d > 1) {
        M <- 2 * (Diag(D) %*% A %*% t(u * dlogft) %*% u.freq + 
            Diag(D) %*% A %*% t(u * dlogT. * dt.dQ) %*% u.freq)
        DA <- M[!lower.tri(M, diag = TRUE)]
    }
    else DA <- NULL
    M <- (A %*% t(u * dlogft) %*% u.freq %*% tA + A %*% t(u * 
        dlogT. * dt.dQ) %*% u.freq %*% tA)
    if (d > 1) 
        DD <- diag(M) + 0.5 * n/D
    else DD <- as.vector(M + 0.5 * n/D)
    grad <- (-2) * c(Dbeta[stack.index[,-(d+1)]], DD * (-2 * D), DA)
    if (is.na(fixed.df)) {
        dlogft.ddf <- 0.5 * (digamma((df + d)/2) - digamma(df/2) - 
            d/df + (df + d) * Q/((1 + Q/df) * df^2) - log(1 + 
            Q/df))
        eps <- 1e-04
        T.eps <- pt(L * sqrt((df + eps + d)/(Q + df + eps)), 
            df + eps + d)
        dlogT.ddf <- (log(T.eps) - log(T.))/eps
        Ddf <- sum((dlogft.ddf + dlogT.ddf) * freq)
        grad <- c(grad, -2 * Ddf * df)
    }
    if (trace) 
        cat("mt.dev.grad: norm is ", sqrt(sum(grad^2)), "\n")
    return(grad)
}#end of mt.dev.grad
#
#  multinomRob
#
#  Walter R. Mebane, Jr.
#  Cornell University
#  http://macht.arts.cornell.edu/wrm1/
#  wrm1@macht.arts.cornell.edu
#
#  Jasjeet Singh Sekhon 
#  Harvard University
#  http://jsekhon.fas.harvard.edu/
#  jsekhon@fas.harvard.edu
#
#  $Id: multinomTanh.R,v 1.10 2004/02/19 02:13:11 wrm1 Exp $
#
#
# mGNtanh, multinomTanh and robustified.leverage
#
#

## mGNtanh:  Gauss-Newton tanh estimator, for overdispersed grouped multinomial GLM
##   Y:  matrix of (overdispersed and contaminated) multinomial counts
##   Ypos:  matrix indicating which in Y are counts (TRUE) and which are not (FALSE).
##   Xarray:  array of regressors,
##      dim(Xarray) = c(n observations, n parameters, n categories)
##   xvec:  vector to indicate all the coefficient parameters in the model
##      (parms by ncats):
##      It has a 1 for an estimated parameter, an integer >1 for an estimated
##      parameter constrained equal to another estimated parameter (all
##      parameters constrained to be equal to one another have the same integer
##      value in xvec) and a 0 otherwize.
##      example:
##      > xvec
##           [,1] [,2] [,3] [,4] [,5]
##      [1,]    1    1    1    1    0
##      [2,]    1    1    1    1    0
##      [3,]    1    1    1    1    0
##      [4,]    1    1    1    1    0
##   tvec: parms by ncats matrix
##   jacstack:  array of regressors,
##      dim(jacstack) = c(n observations, n UNIQUE parameters, n categories)
mGNtanh <- function(bstart, sigma2, resstart,
                      Y, Ypos, Xarray, xvec, tvec,
                      jacstack,itmax=100,print.level=0) {

  ## mcholeskyL:  lower-tri Cholesky matrix for multinomial;
  ##   p is vector of probs,
  mcholeskyL <- function(p) {
    n <- length(p);
    q <- rep(0,n);
    L <- diag(rep(1,n));
    for (i in 1:n) {
      if (i>1) L[i,1:(i-1)] <- -p[i]/q[1:(i-1)];
      if (i<n) q[i] <- 1 - sum(p[1:i]);
    }
    return(L);
  }

  ## mcholeskyLinv:  lower-tri inverse Cholesky matrix for multinomial;
  ##   p is vector of probs,
  mcholeskyLinv <- function(p) {
    n <- length(p);
    q <- rep(0,n);
    Linv <- diag(rep(1,n));
    for (i in 1:n) {
      if (i>1) Linv[i,1:(i-1)] <- p[i]/q[i-1];
      if (i<n) q[i] <- 1 - sum(p[1:i]);
    }
    return(Linv);
  }

  ## mcholeskyD:  Cholesky diagonal;
  ##   p is vector of probs,
  mcholeskyD <- function(p) {
    n <- length(p);
    q <- d <- rep(0,n);
    for (i in 1:n) q[i] <- 1 - sum(p[1:i]);
    d[1] <- p[1]*q[1];
    if (n>2) d[2:n] <- p[2:n]*q[2:n]/q[1:(n-1)];
    return(d);
  }

  ## probfunc: matrix of estimated probabilities
  probfunc <-
    function(Y, Ypos, Xarray, tvec) {
      nobs <- dim(Y)[1]
      ncats <- dim(Y)[2]
      eta <- matrix(0,nobs,ncats)
      for (j in 1:ncats) {
        useobs <- Ypos[,j];
        eta[useobs,j] <- exp(Xarray[useobs,,j] %*% tvec[,j]);
      }
      return( c(1/(eta %*% rep(1,ncats))) * eta )
    }
  ## scorefunc:  score matrix
  scorefunc <-
    function(Ypos, nobs, nparms, phat, N, presmat, wmat, jacstack) {
      scoremat <- matrix(0,nparms,nobs);
      for (i in 1:nobs) {
        usecats <- Ypos[i,];
        nlesscats <- sum(usecats);
        pc <- mcholeskyL(phat[i,usecats]) ;
        pci <- mcholeskyLinv(phat[i,usecats]) ;
        adj <- t(pci) %*% diag(c(wmat[i,1:(nlesscats-1)],1)) %*% t(pc)
        ##     adj <- t(pci) %*% diag(c(wmat[i,1:(nlesscats-1)],1))
        scoremat[,i] <- N[i] * presmat[i,usecats] %*% adj %*% t(jacstack[i,,usecats]) ;  ## weighted
        ##     scoremat[,i] <- N[i] * presmat[i,usecats] %*% t(jacstack[i,,usecats]) ;  ## unweighted
      }
      return( scoremat )
    }
  ## hessianfunc:  hessian matrix with simple weights
  hessianfunc <-
    function(Ypos, nobs, ncats, nparms, phat, N, wmat, jacstack) {
      H <- matrix(0,nparms,nparms)
      for (i in 1:nobs) {
        usecats <- Ypos[i,];
        nlesscats <- sum(usecats);
        pc <- mcholeskyL(phat[i,usecats]) ;
        pD <- mcholeskyD(phat[i,usecats]) ;
        wpvmat <- pc %*% diag(pD * c(wmat[i,1:(nlesscats-1)],1)) %*% t(pc) ;  ## matrix-weighted varmat
        ##     wpvmat <- diag(pD * c(wmat[i,1:(nlesscats-1)],1)) ;  ## matrix-weighted varmat
        ##     wpvmat <- diag(phat[i,usecats])-outer(phat[i,usecats],phat[i,usecats]);  ## unweighted
        H0 <- N[i] * wpvmat;
        H <- H + jacstack[i,,usecats] %*% H0 %*% t(jacstack[i,,usecats])
      }
      return( H )
    }
  ## hessianfunc2:  mean hessian matrix with weights squared
  hessianfunc2 <-
    function(Ypos, nobs, ncats, nparms, phat, N, wmat, jacstack) {
      H <- matrix(0,nparms,nparms)
      for (i in 1:nobs) {
        usecats <- Ypos[i,];
        nlesscats <- sum(usecats);
        pc <- mcholeskyL(phat[i,usecats]) ;
        pD <- mcholeskyD(phat[i,usecats]) ;
        wpvmat <- pc %*% diag(pD * c(wmat[i,1:(nlesscats-1)]^2,1)) %*% t(pc) ;  ## matrix-weighted varmat
        ##     wpvmat <- diag(pD * c(wmat[i,1:(nlesscats-1)]^2,1)) ;  ## matrix-weighted varmat
        ##     wpvmat <- diag(phat[i,usecats])-outer(phat[i,usecats],phat[i,usecats]);  ## unweighted
        H0 <- N[i] * wpvmat;
        H <- H + jacstack[i,,usecats] %*% H0 %*% t(jacstack[i,,usecats])
      }
      return( H / nobs)
    }
  ## resfunc:  orthogonalized and standardized (for multinomial covariance) resids
  resfunc <-
    function(Y, Ypos, Xarray, tvec) {
      if (all(Ypos)) {
        r <- res.std(Y, c(Y %*% rep(1,dim(Y)[2])), probfunc(Y, Ypos, Xarray, tvec));
      }
      else {
        nobs <- dim(Y)[1];
        ncats <- dim(Y)[2];
        r <- matrix(0, nobs, ncats-1);
        phat <- probfunc(Y, Ypos, Xarray, tvec);
        hasall <- apply(Ypos, 1, sum) == ncats;
        nobsall <- sum(hasall);
        if (nobsall > 0) {
          Yuse <- matrix(Y[hasall,], nobsall, ncats);  # in case nobsall == 1
          puse <- matrix(phat[hasall,], nobsall, ncats);
          r[hasall,] <- res.std(Yuse, c(Yuse %*% rep(1,ncats)), puse);
        }
        hasless <- (1:nobs)[!hasall];
        for (i in hasless) {  # orthostd resids go into r[i,1:(nlesscats-1)]
          usecats <- Ypos[i,];
          nlesscats <- sum(usecats);
          Yuse <- matrix(Y[i,usecats], 1, nlesscats);
          puse <- matrix(phat[i,usecats], 1, nlesscats);
          r[i,1:(nlesscats-1)] <- res.std(Yuse, c(Yuse %*% rep(1,nlesscats)), puse);
        }
      }
      return( r );
    }
  psifunc <-
    function(arg) {
      ## Hampel, Rousseeuw and Ronchetti 1981.  constants are from Table 2, p. 645
      ##                                                                          effic.
      ## c <- 3.0;    k <- 5.0;    A <- 0.680593;    B <- 0.769313;    d <- 1.470089;  ## 87%
      c <- 4.0;    k <- 5.0;    A <- 0.857044;    B <- 0.911135;    d <- 1.803134;  ## 97%
      return(
             ifelse(abs(arg)<d, arg,
                    ifelse(abs(arg)<c,
                           sqrt(A*(k-1))*tanh(sqrt((k-1)*B^2/A)*(c-abs(arg))/2)*sign(arg), 0))
             )
    }
  ## compute weight matrix
  weights <-
    function(Y, Ypos, Xarray, tvec, sigma2,ncats) {
      res <- resfunc(Y,Ypos,Xarray,tvec);
      sres <- res/sqrt(sigma2);
      ipsi <- psifunc(sres);
      wNA <- w <- ifelse(sres==0,1,ipsi/sres);
      if (!all(Ypos)) {  # set weights for nonexistent categories to zero
        npos <- apply(Ypos,1,sum);
        nobs <- dim(Y)[1];
        ncats <- dim(Y)[2];
        nwcats <- dim(w)[2];
        for (i in 1:nobs) {
          if (npos[i] < ncats) {
            w[i,npos[i]:nwcats] <- 0;
            wNA[i,npos[i]:nwcats] <- NA;
          }
        }
      }
      list(w=w,ipsi=ipsi,wNA=wNA)
    }
  ## check convergence
  converged <-
    function(bnew,bold) {
      return( sqrt(sum((bnew-bold)^2)) < 1e-6*(sqrt(sum(bold^2)) + 1e-4) )
    }
  converged2 <-
    function(snew,sold) {
      return( sqrt(sum((snew-sold)^2)) <= 1e-8*sqrt(sum(sold^2)) )
    }

  ## begin data computations
  Y[!Ypos] <- 0;  # ensure noncounts are set to zero, for convenience
  ncats <- dim(Y)[2]
  sres <- resstart/sqrt(sigma2);
  wmat <- ifelse(sres==0,1,psifunc(sres)/sres);
  bvec <- bstart;
  ## begin definition of variables used in GNstep that do not change over iterations
  nobs <- dim(Y)[1]
  ncats <- dim(Y)[2]
  catidx <- 1:(ncats-1);
  tvars.total <- dim(Xarray)[2]
  tvunique <- dim(jacstack)[2]
  mvec <- c(Y %*% rep(1,dim(Y)[2]));
  propmat <- Y / mvec;  ## transform observed counts to proportions
  ## end definition of variables used in GNstep that do not change over iterations

  LogLik <-
    function(Y,Ypos,wmatS,ipmatS,mvecS) {
      LL <- 0;  LLu <- 0;
      for (i in 1:nobs) {
        usecats <- Ypos[i,];
        nlesscats <- sum(usecats);
        pc <- mcholeskyL(ipmatS[i,usecats]) ;
        pD <- mcholeskyD(ipmatS[i,usecats]) ;
        dnum <- diag( pc %*% diag(pD * c(wmat[i,1:(nlesscats-1)],0)) %*% t(pc) );
        pvec <- ipmatS[i,usecats];
        dden <- diag( diag(pvec)-outer(pvec,pvec) );
        adj <- dnum/dden;
        LL <- LL - sum(adj * Y[i,usecats] * log(pvec));  ## negative loglikelihood
        LLu <- LLu - sum(Y[i,usecats] * log(pvec));  ## negative loglikelihood
      }
      ##   print(paste("unweighted:",LLu,"; weighted:",LL))
      return(list(LL=LL,LLu=LLu));
    }

  ## Newton algorithm given data and a weight matrix (wmat)
  GNstep <-
    function(bvec,wmatGN,sigma2GN,itmaxGN=100,print.level=0) {
      tvec <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,bvec, ncats,tvars.total);
      itersGN <- 0;
      ## patterned after the Gauss-Newton algorithm in Gallant 1987, 28-29
      for (iGN in 1:itmaxGN) {
        itersGN <- itersGN + 1;
        bprev2 <- bvec;
        
        ipmat <- probfunc(Y, Ypos, Xarray, tvec);
        presmat <- propmat - ipmat;
        loglik <- LogLik(Y,Ypos,wmatGN,ipmat,mvec)$LL;
        if (iGN==1 & (print.level > 1) )
          print(paste("mGNtanh: -loglik initial:",loglik));
        score <-
          scorefunc(Ypos, nobs, tvunique, ipmat, mvec, presmat, wmatGN, jacstack);
        hess2 <-
          hessianfunc2(Ypos, nobs, ncats, tvunique, ipmat, mvec, wmatGN, jacstack);
        posdef <- all(eigen(hess2, symmetric=TRUE, only.values=TRUE)$values > 0);
        gradient <- (score %*% rep(1,nobs)) / sqrt(sigma2GN) / nobs ;
        if (!posdef) {
          convflag <- FALSE;
          break;  ## quit if Hessian is not positive definite
        }
        bdiff <- c(solve(hess2, tol=.Machine$double.eps, LINPACK=TRUE) %*% gradient) ;  ## one Newton step
        ## print(bdiff);
        for (lambda in c(10:6)/10) {
          blambda <- bvec + lambda * bdiff;
          tlambda <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,blambda, ncats,tvars.total);
          plambda <- probfunc(Y, Ypos, Xarray, tlambda);
          logliklambda <- LogLik(Y,Ypos,wmatGN,plambda,mvec)$LL;
          if (!is.na(logliklambda) && logliklambda < loglik) break;
        }
        if (is.na(logliklambda) || logliklambda > loglik) for (lambda in 2^(-(1:45))) {
          blambda <- bvec + lambda * bdiff;
          tlambda <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,blambda, ncats,tvars.total);
          plambda <- probfunc(Y, Ypos, Xarray, tlambda);
          logliklambda <- LogLik(Y,Ypos,wmatGN,plambda,mvec)$LL;
          if (!is.na(logliklambda) && logliklambda < loglik) break;
        }
        if (logliklambda < loglik) bvec <- blambda;
        tvec <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,bvec, ncats,tvars.total);
        if (print.level > 1)  {
          cat("mGNtanh: ibvec: "); print(bvec)
        }
        convflag <- converged(bvec,bprev2) & converged2(logliklambda,loglik);
        ##     convflag <- convflag & all(abs(gradient) < 1e-9);
        if (convflag) break;
      }
      if (!posdef & print.level >= 0) {
        print("mGNtanh: Hessian is not positive definite");
      }
      if (print.level > 1 & posdef) {
        print(paste("mGNtanh: -loglik final: ",logliklambda));
      }
      LL2 <- ifelse(posdef, LogLik(Y,Ypos,wmatGN,plambda,mvec), NA);
      if (print.level > 1 & posdef) {
        print(paste("mGNtanh: -loglik:  weighted,",LL2$LL,";  unweighted,",LL2$LLu));
        print("mGNtanh: gradient:");  print(c(gradient));
        print("mGNtanh: bvec:");  print(bvec);
      }
      information <-
        hessianfunc(Ypos, nobs, ncats, tvunique, ipmat, mvec, wmatGN, jacstack);
      if (all(eigen(information, symmetric=TRUE, only.values=TRUE)$values > 0)) {
        formation <- solve(information, tol=.Machine$double.eps, LINPACK=TRUE);
      }
      else {
        formation <- NA;
      }
      return(
             list(coefficients=bvec, tvec=tvec, formation=formation, score=score,
                  LLvals=LL2, convflag=convflag, iters=itersGN, posdef=posdef) );
    }
  error <- 0;
  iters <- 0;
  for (i in 1:itmax) {
    iters <- iters + 1;
    bprev <- bvec;
    wprev <- wmat;

    tvec <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,bvec, ncats,tvars.total);
    wres <- resfunc(Y,Ypos,Xarray,tvec)*wmat;
    wobs  <- sum(wmat);
    tanhsigma2  <- sum(wres^2)/(wobs-tvunique);

    ## grouped multinomial:  estimate using Newton algorithm
    GNlist <- GNstep(bvec,wmat,tanhsigma2,print.level);
    error <- ifelse(GNlist$posdef,0,32);  ## error == 32 if hessian not posdefinite
    if (print.level > 1)
      print(paste("mGNtanh: number of Newton iterations", GNlist$iters));

    bvec <- GNlist$coeff;
    wlist <- weights(Y,Ypos,Xarray,GNlist$tvec,sigma2,ncats);
    wmat <- wlist$w;
    wmatNA <- wlist$wNA;
    if (converged(bvec,bprev) & converged(wmat,wprev)) break;
  }

  opg <- GNlist$score %*% t(GNlist$score) ;
  obsformation <- GNlist$formation ;
  rcovmat <- obsformation %*% opg %*% obsformation;

  if (print.level > 1) {  
    if (length(obsformation)==1 && obsformation==NA) {
      print(paste("mGNtanh: hessian determinant:",NA));
    }
    else {
      print(paste("mGNtanh: hessian determinant:",
        det(solve(obsformation, tol=.Machine$double.eps, LINPACK=TRUE))));
    }
    print(paste("mGNtanh: OPG determinant:", det(opg)));
    print(paste("mGNtanh: tanh sigma^2:", tanhsigma2));
  }

#  if (tanhsigma2 > sigma2) error <- error + 1;  ## error: tanh sigma2 > LQD sigma2
  if (sum(wmat) < nobs*(ncats-1)/2) error <- error + 2;  ## error: wgts are too small

  ## table of returned error values (indicated values add to give total error)
  ## 0    no errors
  ## 1   tanh sigma2 > LQD sigma2
  ## 2   sum of weights < nobs*(ncats-1)/2
  ## 32  Hessian not positive definite in the final Newton step

  return(
         list(coefficients=GNlist$tvec, coeffvec=GNlist$coeff, dispersion=sigma2,
              w=wlist$wNA, psi=wlist$ipsi,
              A=opg/tanhsigma2, B=obsformation, covmat=rcovmat,
              iters=iters, error=error,
              GNlist=GNlist, tanhsigma2=tanhsigma2,
              Y=Y, Ypos=Ypos, probmat=probfunc(Y, Ypos, Xarray, GNlist$tvec),
              jacstack=jacstack, Xarray=Xarray)
         )
}

multinomTanh <- function (Y, Ypos, X, jacstack, xvec, tvec, pop, s2,
                     xvar.labels, choice.labels,
                     print.level=0) {

    nobs  <- dim(Y)[1];
    ncats <- dim(Y)[2];
    nvars <- dim(X)[2];
    nvars.unique <- sum(xvec == 1) + length(unique(xvec[xvec>1]));


    beta.vector <- vector(mode="numeric", length=nvars.unique);
    beta.vector <- mnl.xvec.mapping(forward=TRUE,xvec, tvec, beta.vector,
                                   ncats,nvars);

    residuals <- residual.generator(tvec,Y,Ypos,X,pop);

    mtanh <-
      mGNtanh(beta.vector, s2, residuals$Sres,
               Y, Ypos, X, xvec, as.matrix(tvec),
               jacstack,print.level);

    tvec <- mnl.xvec.mapping(forward=FALSE,xvec,tvec,mtanh$coeffvec,
                             ncats,nvars);

    se.vec <- sqrt(diag(mtanh$covmat));
    se  <- xvec;
    se  <- as.data.frame(mnl.xvec.mapping(forward=FALSE,xvec,se,se.vec,
                                          ncats,nvars));

    se.opg.vec  <- sqrt(diag(ginv(mtanh$A)));
    se.opg  <- xvec;
    se.opg  <- as.data.frame(mnl.xvec.mapping(forward=FALSE,xvec,se.opg,se.opg.vec,
                                              ncats,nvars));

    se.hes.vec <- sqrt(diag(mtanh$B));
    se.hes  <- xvec;
    se.hes  <- as.data.frame(mnl.xvec.mapping(forward=FALSE,xvec,se.hes,se.hes.vec,
                                              ncats,nvars));        

    row.names(se) <- xvar.labels;
    names(se)     <- choice.labels;
    mtanh$se      <- se;

    row.names(se.opg) <- xvar.labels;
    names(se.opg)     <- choice.labels;
    mtanh$se.opg      <- se.opg;

    row.names(se.hes) <- xvar.labels;
    names(se.hes)     <- choice.labels;
    mtanh$se.hes      <- se.hes;

    Hdiag <- robustified.leverage(tvec, Y, Ypos, X, pop, ifelse(mtanh$w>0,1,0),jacstack);
    w.Hdiag <- as.data.frame(matrix(c(as.vector(1:nobs),
                                     signif(mtanh$w), signif(Hdiag)),ncol=(ncats-1)+(ncats-1)+1));
    names(w.Hdiag) <-
      c("name",paste("weights:",choice.labels[1:ncats-1],sep=""),
        paste("Hdiag:",choice.labels[1:ncats-1],sep=""));
    #cat("mtanh: weights, Hdiag (by choices)\n");

    sigma2 <- mtanh$disp;
    sigma  <- sqrt(sigma2);    

    cr <- fn.region.results(tvec, Y, Ypos, X, pop, sigma, Hdiag);

    mtanh$coef <- tvec;
    weights <- w.Hdiag[,c(1:ncats)];
    j  <- ncats
    Hdiag   <- w.Hdiag[,c(1,(j+1):(j+j-1))];

    return( list(mtanh= mtanh,
                 weights=weights,
                 Hdiag=Hdiag,
                 cr   = cr,
                 tvec =tvec));
  } #multinomTanh


robustified.leverage <- function (tvec, Y, Ypos, Xarray, m, Win,jacstack) {
#tvec <- mout$mtanh$coef
#Xarray <- X;
#m <- TotalY;
#W <- ifelse(mtanh$w>0,1,0);
#Win  <- mtanh$w

  obs  <- dim(jacstack)[1];
  W  <- cbind(Win,rep(1,obs));
  
  nobs  <- dim(Y)[1];
  nvars <- dim(Xarray)[2];
  ncats <- dim(Xarray)[3];

  y.prob <- mnl.probfunc(Y, Ypos, Xarray, tvec);
  allYpos <- all(Ypos);
  if (!allYpos) {  # put probs for existing categories in y.prob[i,1:npos]
    npos <- apply(Ypos,1,sum);
    for (i in 1:nobs) {
      if (npos[i] < ncats) {
        y.prob[i,1:npos[i]] <- y.prob[i,Ypos[i,]];
        y.prob[i,(npos[i]+1):ncats] <- 0;
      }
    }
  }
  
  Hdiag <- matrix(0,nobs,ncats-1);
  summat <- matrix(0,ncats,ncats)
  for (i in 1:(ncats-1)) summat[(i+1):ncats,i] <- 1;
  p <- y.prob;
  q <- p %*% summat;

  #tanabe sagae '92
  L.gen <- function(p,q,ncats) {
    L <- matrix(0,ncats,ncats)
    for (j1 in 1:ncats) {
      for (j2 in 1:ncats) {
        if (j1 > j2)  { L[j1,j2]  <- -p[j1]/q[j2]; }
        if (j1 == j2) { L[j1,j2] <- 1; }
        if (j1 < j2)  { L[j1,j2] <- 0; }
      }
    }
    return(L)
  } #end of L.gen
  
  invL <- function(p,q,ncats) {
    Linv <- matrix(0,ncats,ncats)
    for (j1 in 1:ncats) {
      for (j2 in 1:ncats) {
        if (j1 > j2)  { Linv[j1,j2]  <- p[j1]/q[j1-1]; }
        if (j1 == j2) { Linv[j1,j2] <- 1; }
        if (j1 < j2)  { Linv[j1,j2] <- 0; }
      }
    }
    return(Linv)
  } #end of invL

# d:  nonzero values in the diagonal matrix of the Cholesky decomposition
  # note:  d[j,i]==0 if fewer than i+1 alternatives exist for obs j
  d <- matrix(0,obs,ncats-1)
  for (i in 1:(ncats-1)) {
    if (i==1) d[,i] <- p[,i]*q[,i];
    if (i>1)  d[,i] <- p[,i]*q[,i]/q[,i-1];
  }  

  Center  <- matrix(0,dim(jacstack)[2],dim(jacstack)[2]);
  for (i in 1:obs)
    {
      V <- matrix(0,nrow=ncats,ncol=ncats);
      V <- diag(c( ifelse(d[i,]>0, 1/sqrt(m[i]*d[i,]), 0) ,0));
      #USE THIS??? V <- V*ncats^2;

      L  <-  L.gen(p[i,],q[i,],ncats);
      
      if (!allYpos) {  # put jacstack for existing categories in smp[,1:npos]
        smp <- jacstack[i,,];
        if (npos[i] < ncats) {
          smp[,1:npos[i]] <- smp[,Ypos[i,]];
          smp[,(npos[i]+1):ncats] <- 0;
        }
        C0  <- smp %*% t(L) %*% V %*% diag(W[i,]) %*% V %*% t(L) %*% t(smp);
      }
      else {
        C0  <- jacstack[i,,] %*% t(L) %*% V %*% diag(W[i,]) %*%
                 V %*% t(L) %*% t(jacstack[i,,]);
      }
      Center  <- Center + C0;
    } #end of i

  invCenter  <- ginv(Center);

  for (i in 1:obs)
    {
      V <- matrix(0,nrow=ncats,ncol=ncats);
      V <- diag(c( ifelse(d[i,]>0, 1/sqrt(m[i]*d[i,]), 0) ,0));
      L  <-  L.gen(p[i,],q[i,],ncats);

      if (!allYpos) {  # put jacstack for existing categories in smp[,1:npos]
        smp <- jacstack[i,,];
        if (npos[i] < ncats) {
          smp[,1:npos[i]] <- smp[,Ypos[i,]];
          smp[,(npos[i]+1):ncats] <- 0;
        }
        Hdiag[i,]  <- diag(V %*% t(L) %*% t(smp) %*%
                        invCenter %*% smp %*% t(L) %*% V)[1:(ncats-1)];
      }
      else {
        Hdiag[i,]  <-
          diag(V %*% t(L) %*% t(jacstack[i,,]) %*%
            invCenter %*% jacstack[i,,] %*% t(L) %*% V)[1:(ncats-1)];
      }
    } #end of i
  
# put the negative forecasting variance adjustment values in Hdiag[w==0]
  for (j in 1:(ncats-1)) {
    sindx  <- Win[,j]==0;
    Hdiag[sindx,j] <- -Hdiag[sindx,j];
  }
  
  return(Hdiag);
}    #end of robustified.leverage

#calculate region results (orthogonal)
fn.region.results <- function (tmp.vec, Y, Ypos, X, TotalY, sigma, Hdiag) {
  y.prob      <- mnl.probfunc(Y,Ypos,X,tmp.vec);

  if (all(Ypos)) {
    Sres.raw <- res.std(Y, TotalY, y.prob);
  }
  else {
    nobs <- dim(Y)[1];
    ncats <- dim(Y)[2];
    Sres.raw <- matrix(0, nobs, ncats-1);
    hasall <- apply(Ypos, 1, sum) == ncats;
    nobsall <- sum(hasall);
    if (nobsall > 0) {
      Yuse <- matrix(Y[hasall,], nobsall, ncats);  # in case nobsall == 1
      puse <- matrix(y.prob[hasall,], nobsall, ncats);
      Sres.raw[hasall,] <- res.std(Yuse, TotalY[hasall], puse);
    }
    hasless <- (1:nobs)[!hasall];
    for (i in hasless) {
      usecats <- Ypos[i,];
      nlesscats <- sum(usecats);
      ocats <- 1:(nlesscats-1);
      Yuse <- matrix(Y[i,usecats], 1, nlesscats);
      puse <- matrix(y.prob[i,usecats], 1, nlesscats);
      Sres.raw[i,ocats] <- res.std(Yuse, TotalY[i], puse);
    }
  }

  Hmax <- .9;

  standard <- Sres.raw / sigma;
  student  <- Sres.raw / (sigma * sqrt(1-ifelse(Hdiag<Hmax,Hdiag,Hmax)));

  return(list(pred=y.prob,student=student,standard=standard));
} #end robustified leverage
#
#  multinomRob
#
#  Walter R. Mebane, Jr.
#  Cornell University
#  http://macht.arts.cornell.edu/wrm1/
#  wrm1@macht.arts.cornell.edu
#
#  Jasjeet Singh Sekhon 
#  Harvard University
#  http://jsekhon.fas.harvard.edu/
#  jsekhon@fas.harvard.edu
#
#  $Id: newtanh8a.R,v 1.10 2004/02/17 05:17:01 wrm1 Exp $
#

## *** NOTE:  mGNtanh definition moved into multinomTanh.R ***
#
#  multinomRob
#
#  Walter R. Mebane, Jr.
#  Cornell University
#  http://macht.arts.cornell.edu/wrm1/
#  wrm1@macht.arts.cornell.edu
#
#  Jasjeet Singh Sekhon 
#  Harvard University
#  http://jsekhon.fas.harvard.edu/
#  jsekhon@fas.harvard.edu
#
#  $Id: resstd2.R,v 1.5 2004/03/04 02:08:17 wrm1 Exp $
#

# probfunc: matrix of estimated probabilities
mnl.probfunc <-  function(Y, Ypos, Xarray, tvec) {
  nobs <- dim(Y)[1]
  ncats <- dim(Y)[2]
  eta <- matrix(0,nobs,ncats)
  for (j in 1:ncats) {
    useobs <- Ypos[,j];
    eta[useobs,j] <- exp(Xarray[useobs,,j] %*% tvec[,j]);
  }
  return( c(1/(eta %*% rep(1,ncats))) * eta )
}#end of mnl.probfunc

## residual.generator:  raw and median-centered ortho-standardized residuals
residual.generator <- function (tvec, Y, Ypos, X, m)
  {

    y.prob <- mnl.probfunc(Y, Ypos, X, tvec);

    if (all(Ypos)) {
      Sres.raw <- res.std(Y, m, y.prob);
      Sres <- Sres.raw - median(Sres.raw);
    }
    else {
      nobs <- dim(Y)[1];
      ncats <- dim(Y)[2];
      Sres.raw <- matrix(0, nobs, ncats-1);
      hasall <- apply(Ypos, 1, sum) == ncats;
      nobsall <- sum(hasall);
      if (nobsall > 0) {
        Yuse <- matrix(Y[hasall,], nobsall, ncats);  # in case nobsall == 1
        puse <- matrix(y.prob[hasall,], nobsall, ncats);
        Sres.raw[hasall,] <- res.std(Yuse, c(Yuse %*% rep(1,ncats)), puse);
      }
      hasless <- (1:nobs)[!hasall];
      Sres.use <- matrix(TRUE, nobs, ncats-1);
      for (i in hasless) {
        usecats <- Ypos[i,];
        nlesscats <- sum(usecats);
        ocats <- 1:(nlesscats-1);
        Yuse <- matrix(Y[i,usecats], 1, nlesscats);
        puse <- matrix(y.prob[i,usecats], 1, nlesscats);
        Sres.raw[i,ocats] <- res.std(Yuse, c(Yuse %*% rep(1,nlesscats)), puse);
        Sres.use[i,nlesscats:(ncats-1)] <- FALSE;
      }
      Sres <- Sres.raw - median(Sres.raw[Sres.use]);
      Sres[!Sres.use] <- 0;
    }

    return(list(Sres.raw=Sres.raw,Sres=Sres));
  }

## res.std: compute orthogonalized and standardized residuals
##  based on Tanabe and Sagae (1992 JRSSB)
res.std <- function(y,m,p, print.level=0)
  {
    obs    <- dim(y)[1]
    ncats  <- dim(y)[2]
    summat <- matrix(0,ncats,ncats)
    for (i in 1:(ncats-1)) summat[(i+1):ncats,i] <- 1;
    q <- p %*% summat
    ## d:  nonzero values in the diagonal matrix of the Cholesky decomposition
    d <- matrix(0,obs,ncats-1)
    for (i in 1:(ncats-1)) {
      if (i==1) d[,i] <- p[,i]*q[,i];
      if (i>1)  d[,i] <- p[,i]*q[,i]/q[,i-1];
    }
    ## r: raw residuals
    r <- y-m*p;
    summat <- matrix(0,ncats,ncats)
    for (i in 1:ncats) summat[i,i:ncats] <- 1;
    rsum <- r %*% summat
    ## Or:  orthogonalized residuals
    Or <- matrix(0,obs,ncats-1)
    for (i in 1:(ncats-1)) {
      if (i==1) Or[,i] <- r[,i];
      if (i>1)  Or[,i] <- r[,i] + rsum[,i-1]*p[,i]/q[,i-1];
    }
    ## Sr: standardized residuals
    Sr <- Or/sqrt(m*d);

    if (print.level > 0) {
      cat("res.std: Or\n");
      print(Or);
      cat("res.std: m\n");
      print(m);
      cat("res.std: d\n");
      print(d); 
    }
    return(Sr);
  }

ResStd  <- function(nobs, ncats, nvars.total, nvars.unique,
                    tvec, Y, X, weights)
  {
    SresRaw  <- matrix(0,nrow=obs,ncol=(ncats-1));
    SresRaw  <- .Call("InResStd", as.integer(nobs),  as.integer(ncats),
                      as.integer(nvars.total), as.integer(nvars.unique),
                      as.real(tvec), as.real(Y), as.real(X),
                      as.real(weights), as.real(SresRaw),
                      PACKAGE="multinomRob");
    return(SresRaw);
  }

kth.smallest  <- function(SortVector, obs, k)
  {
    if (obs < k) {
      print("ERROR! obs < k in kth.smallest\n");
      return(-1);
    }
    return(.Call("kthSmallest",
                 as.real(SortVector), as.integer(obs), as.integer(k),
                 PACKAGE="multinomRob"));    
  } #end of kth.smallest

#
#  multinomRob
#
#  Walter R. Mebane, Jr.
#  Cornell University
#  http://macht.arts.cornell.edu/wrm1/
#  wrm1@macht.arts.cornell.edu
#
#  Jasjeet Singh Sekhon 
#  Harvard University
#  http://jsekhon.fas.harvard.edu/
#  jsekhon@fas.harvard.edu
#
#  $Id: simulation.functions2.R,v 1.2 2004/02/14 06:28:13 wrm1 Exp $
#

# Generate single random Multinomial(n,pr)
rmultinomial<-function(n=5, pr=c(0.5,0.5), long=FALSE) {
  k<-length(pr)
  if (abs(1-sum(pr))>0.000001)
   stop("(rmultinomial): parameter pr must be the k probabilities (summing to 1)")

  if(long) {
    y<-runif(n, 0, 1)
    p<-cumsum(pr)
    Seq<-1:n
    x<-sapply(y, function(y, Seq, p) {Seq[y <= p][1]}, Seq=Seq, p=p)
  } else {
    x<-rep(NA,k)
    p<-pr/c(1,(1-cumsum(pr[1:(k-1)])))
    for (i in 1:(k-1)) {
      if (n==0) {
        x[i]<-0
        if (i==k-1) x[k]<-0
        next
      }
      y<-rbinom(1,n,p[i])
      x[i]<-y
      if (i==k-1) x[k]<-n-y
      n<-n-y
    }
  }
  return(x)
}

# Generate n random Dirichlet's
rdirich<-function(n,alphavec) {
  # Gelman, et al., Bayesian Data Analysis, Chapman & Hall, 1995, p.482
  k<-length(alphavec)
  p<-matrix(rgamma(n*k,alphavec),n,k,byrow=T)
  sm<-matrix(apply(p,1,sum),n,k)
  return(p/sm)
}


#
#  multinomRob
#
#  Walter R. Mebane, Jr.
#  Cornell University
#  http://macht.arts.cornell.edu/wrm1/
#  wrm1@macht.arts.cornell.edu
#
#  Jasjeet Singh Sekhon 
#  Harvard University
#  http://jsekhon.fas.harvard.edu/
#  jsekhon@fas.harvard.edu
#
#  $Id: spec.R,v 1.4 2004/02/14 22:21:35 wrm1 Exp $
#

# functions to interpret model formula specifications and build data to analyze

# get.xdata:
# Return model matrix corresponding to the formula in formul
# 
get.xdata <- function(formul, datafr) {
  t1 <- terms(formul, data=datafr);
  if (length(attr(t1, "term.labels"))==0 & attr(t1, "intercept")==0) {
    m <- NULL;  # no regressors specified for the model matrix
  }
  else {
    m <- model.matrix(formul, data=datafr);
  }
  return(m);
}

# get.ydata:
# Return response vector corresponding to the formula in formul
# 
get.ydata <- function(formul, datafr) {
  t1 <- terms(formul, data=datafr);
  if (length(attr(t1, "response"))==0) {
    m <- NULL;  # no response variable specified
  }
  else {
    m <- model.response(model.frame(formul, data=datafr));
  }
  return(m);
}

# get.xy:
#  Return response vectors and model matrices corresponding the formulas
#  in formlist, along with the variable names and the number of variables
#  in the model matrix for each response.  Each response is a distinct column
#  in returned matrix Y, and each model matrix is along a dimension of the
#  returned array X.  X contains 0 for columns that are not used for a
#  category of the response (corresponding to a column in Y).
#  Factor variables are expanded in X, as binary {0,1} dummy variables, reduced as
#  appropriate to match other factors and the intercept (e.g., with one factor
#  variable (facvar) and an intercept, there are length(levels(facvar))-1 dummy
#  variables, while with one factor variable and no intercept there are
#  length(levels(facvar)) dummy variables).
#  Missing data (NA) in any variable causes the entire observation to be deleted
#  (listwise deletion).
#  A matrix ypos is created to indicate which elements of the response matrix Y
#  originally had negative values:  ypos <- Y >= 0.  The ypos matrix is used to
#  indicate the number of outcome alternatives for each observation.  If there
#  are fewer than two alternatives for an observation, that observation is deleted.
#  dim(Y) == c(nobs,ncats)
#  dim(X) == c(nobs, max(c(1,xlengths)), ncats)
#  length(xlengths) == ncats
#  length(ynames) == ncats & is.character(ynames[i])
#  length(xnames) == ncats & is.list(xnames) & length(xnames[[i]]) == xlengths[i]
#     & is.character(xnames[[i]][j])
# 
get.xy <- function(formlist, datafr, print.level=0) {
  ncats <- length(formlist);
  nobs <- dim(datafr)[1];  # assume all variables are the same length
  Y <- NULL;
  Xlist <- list();
  ynames <- rep("", ncats);
  xlengths <- rep(0, ncats);
  # find negative response values and missing response or regressor variable values
  # if y[j,i] is negative, ignore missing data in x[j,i]
  ypos <- matrix(FALSE, nobs, ncats);  #  in ypos, FALSE is y<0, TRUE is y>=0 or is.na(y)
  ymiss <- matrix(FALSE, nobs, ncats);  #  in ymiss, TRUE is NA, FALSE is not NA
  xmiss <- matrix(FALSE, nobs, ncats);  #  in xmiss, TRUE is NA, FALSE is not NA
  for (i in 1:ncats) {
    ti <- terms(formlist[[i]], data=datafr);
    ynames[i] <- as.character(attr(ti, "variables")[[1 + attr(ti, "response")]]) ;
    ypos[,i] <- datafr[[ ynames[i] ]] >= 0;
    ypos[,i] <- ifelse(is.na(ypos[,i]), TRUE, ypos[,i]);
    ymiss[,i] <- is.na(datafr[[ ynames[i] ]]);
    xni <- as.character(attr(ti, "term.labels")) ;
    if (length(xni) >= 1) {
      for (ii in 1:length(xni)) {
        xmiss[,i] <- xmiss[,i] | is.na(datafr[[ xni[ii] ]]);
      }
    }
  }
  xmiss <- xmiss & ypos;
  # remove any observation for which less than two alternatives exist
  hastwo <- apply(ypos, 1, sum) >= 2;
  if (any(!hastwo)) {
    xmiss[!hastwo,] <- TRUE;
    if (print.level > 0) {
      print("some observations have fewer than two responses.");
    }
  }
  # create vector for listwise deletion
  if (print.level > 0) {
    if (any(ymiss)) {
      print("there is missing response variable data.");
    }
    if (any(xmiss)) {
      print("there is missing regressor variable data.");
    }
  }
  misslist <- apply(ymiss,1,any) | apply(xmiss,1,any);
  if (print.level > 0) {
    if (any(misslist)) {
      print("implementing listwise deletion for missing data and fewer than two responses.");
      cat("deleting observations:  ", c(1:nobs)[misslist], "\n");
    }
    # print(misslist)
  }
  # begin implementation of listwise deletion
  datafr <- datafr[!misslist,];
  ypos <- ypos[!misslist,];
#  nobs <- sum(!misslist);
  # end implementation of listwise deletion

  XYdata <- get.XYdata(formlist, datafr);
  
  dimnames(XYdata$Y) <- list(NULL, XYdata$ynames);
  dimnames(XYdata$X) <- list(NULL, NULL, XYdata$ynames);
  return(list(Y=XYdata$Y, X=XYdata$X, xnames=XYdata$xnames, ynames=XYdata$ynames,
              xlengths=XYdata$xlengths, ypos=ypos));
}

get.XYdata <- function(formlist, datafr) {
  ncats <- length(formlist);
  nobs <- dim(datafr)[1];  # assume all variables are the same length
  Y <- NULL;
  Xlist <- list();
  xlengths <- rep(0, ncats);
  ynames <- rep("", ncats);
  xnames <- list()
  for (i in 1:ncats) {
    ti <- terms(formlist[[i]], data=datafr);
    ynames[i] <- as.character(attr(ti, "variables")[[1 + attr(ti, "response")]]) ;
    Y <- cbind(Y, get.ydata(formlist[[i]], datafr));
    if (length(attr(ti, "term.labels"))==0 & attr(ti, "intercept")==0) {
      Xlist[[i]] <- 0;
      xlengths[i] <- 0;
      xnames[[i]] <- "";
    }
    else {
      Xlist[[i]] <- get.xdata(formlist[[i]], datafr);
      xlengths[i] <- dim(Xlist[[i]])[2];
      xnames[[i]] <- unlist(dimnames(Xlist[[i]])[2]);
    }
  }
  X <- array(0, dim=c(nobs, max(c(1,xlengths)), ncats));
  for (i in 1:ncats) {
    if (xlengths[i] > 0) {
      X[1:nobs, 1:xlengths[i], i] <- Xlist[[i]];
    }
  }
  return(list(Y=Y, X=X, xnames=xnames, ynames=ynames, xlengths=xlengths));
}

get.xynames <- function(formlist, datafr) {
  ncats <- length(formlist);
  ynames <- rep("", ncats);
  xnames <- list()
  for (i in 1:ncats) {
    ti <- terms(formlist[[i]], data=datafr);
    ynames[i] <- as.character(attr(ti, "variables")[[1 + attr(ti, "response")]]) ;
    if (length(attr(ti, "term.labels"))==0 & attr(ti, "intercept")==0) {
      xnames[[i]] <- "";
    }
    else {
      xnames[[i]] <- unlist(dimnames(get.xdata(formlist[[i]], datafr))[2]);
    }
  }
  return(list(xnames=xnames, ynames=ynames));
}
#
#  multinomRob
#
#  Walter R. Mebane, Jr.
#  Cornell University
#  http://macht.arts.cornell.edu/wrm1/
#  wrm1@macht.arts.cornell.edu
#
#  Jasjeet Singh Sekhon 
#  Harvard University
#  http://jsekhon.fas.harvard.edu/
#  jsekhon@fas.harvard.edu
#
#  $Id: zzz.R,v 1.5 2004/02/19 02:13:11 wrm1 Exp $
#

# use .onLoad instead of .First.lib for use with NAMESPACE and R(>= 1.7.0)
.onLoad <- function(lib, pkg) {
  library.dynam(pkg, pkg, lib)
  require(rgenoud)||
    cat("ERROR: library 'rgenoud' is needed by 'multinomRob' and is missing\n")
  require(MASS)   ||
    cat("ERROR: library 'MASS' is needed by 'multinomRob' and is missing\n")  
  require(mvtnorm)||
    cat("ERROR: library 'mvtnorm' is needed by 'multinomRob' and is missing\n")
}#end of .First

#.First.lib <- function(lib, pkg) {
#  library.dynam(pkg, pkg, lib)
#  require(rgenoud)||
#    cat("ERROR: library 'rgenoud' is needed by 'multinomRob' and is missing\n")
#  require(MASS)   ||
#    cat("ERROR: library 'MASS' is needed by 'multinomRob' and is missing\n")  
#  require(mvtnorm)||
#    cat("ERROR: library 'mvtnorm' is needed by 'multinomRob' and is missing\n")
#}#end of .First

.onUnload <- function(libpath) {
   library.dynam.unload("multinomRob", libpath)
}
