```r
library(MASS)
predict.hier.mnl <- function(model, data, nresp=1000) {
  data.model <-
    model.matrix(update(model$formula, 0 ~ .),
                 data = data)[,-1]
  coef.Sigma <- cov.mlogit(model)
  coef.mu <- m2.hier$coef[1:dim(coef.Sigma)[1]]
  draws <- mvrnorm(n=nresp, coef.mu, coef.Sigma)
  shares <- matrix(NA, nrow=nresp, 
                   ncol=nrow(data))
  for (i in 1:nresp) {
    utility <- data.model %*% draws[i, ]
    share = exp(utility) / sum(exp(utility))
    shares[i, ] <- share
  }
  cbind(colMeans(shares), data)
}
```