目录
多元离散选择模型简介样例程序样例程序的问题问题解决
多元离散选择模型简介
常用的离散选择模型有logit模型及probit模型，其区别就是假设不可观测量的分布不同。logit模型假设不可观测项服从Gumbel (Extreme Value Type I) Distribution 。多元logit模型是，多元离散选择模型的一种，适合于效用最大化时的分布选择。
如果决策者i在J项可供选择方案中选择了第j项，那么其效用模型为
以下是基于R语言的一个程序
#引入所需要的包require(foreign)require(nnet)require(ggplot2)require(reshape2)##读取数据集，这里数据集可以来找csv，mysql等ml <- read.dta("http://www.ats.ucla.edu/stat/data/hsbdemo.dta")head(ml) ## 查看数据的前六行with(ml, table(ses, prog)) ## 以表格的形式统计数目##函数tapply（进行分组统计）with(ml, do.call(rbind, tapply(write, prog,function(x) c(M = mean(x), SD = sd(x)))))## 重新排列因子水平，其中以academic作为基础标准ml$prog2 <- relevel(ml$prog, ref = "academic")ml$prog2##训练模型test <- multinom(prog2 ~ ses + write, data = ml)summary(test)#获取z值z <- summary(test)$coefficients/summary(test)$standard.errorsz#获取p值p <- (1 - pnorm(abs(z), 0, 1))*2p## extract the coefficients from the model and exponentiateexp(coef(test))head(pp <- fitted(test))##测试集上检测dses <- data.frame(ses = c("low", "middle", "high"),write = mean(ml$write))#获取预测的概率predict(test, newdata = dses, "probs")dwrite <- data.frame(ses = rep(c("low", "middle", "high"), each = 41),write = rep(c(30:70), 3))## store the predicted probabilities for each value of ses and writepp.write <- cbind(dwrite, predict(test, newdata = dwrite, type = "probs", se = TRUE))## calculate the mean probabilities within each level of sesby(pp.write[, 3:5], pp.write$ses, colMeans)## melt data set to long for ggplot2lpp <- melt(pp.write, id.vars = c("ses", "write"), value.name = "probability")head(lpp) # view first few rows## plot predicted probabilities across write values for## each level of ses facetted by program typeggplot(lpp, aes(x = write, y = probability, colour = ses)) + geom_line() + facet_grid(variable ~ ., scales="free") 样例程序的问题但是，在使用上面的样例，如下，当X的维数较大时，类别较多时，就会出现问题，其主要问题是迭代时权重的问题，为此，我们需要修改源码。
test <- multinom(trainset$car_id2 ~ ., data = trainset[,-c(1,13,14)])以下是容易出现的问题：
Error in nnet.default(X, Y, w, mask = mask, size = 0, skip = TRUE, softmax = TRUE, : too many (11745) weights 问题解决针对问题too many (11745) weights，这类问题，解决的办法是修改源码。即新写一个方法。如下，我修改的部分是这里面的
MaxNWts=12000 ##将程序中的MaxNWts变大就行了 ##估计参数的函数nnetfunction<-function (formula, data, weights, subset, na.action, contrasts = NULL, Hess = FALSE, summ = 0, censored = FALSE, model = FALSE, ...) { class.ind <- function(cl) { n <- length(cl) x <- matrix(0, n, length(levels(cl))) x[(1L:n) + n * (as.integer(cl) - 1L)] <- 1 dimnames(x) <- list(names(cl), levels(cl)) x } summ2 <- function(X, Y) { X <- as.matrix(X) Y <- as.matrix(Y) n <- nrow(X) p <- ncol(X) q <- ncol(Y) Z <- t(cbind(X, Y)) storage.mode(Z) <- "double" z <- .C(VR_summ2, as.integer(n), as.integer(p), as.integer(q), Z = Z, na = integer(1L)) Za <- t(z$Z[, 1L:z$na, drop = FALSE]) list(X = Za[, 1L:p, drop = FALSE], Y = Za[, p + 1L:q]) } call <- match.call() m <- match.call(expand.dots = FALSE) m$summ <- m$Hess <- m$contrasts <- m$censored <- m$model <- m$... <- NULL m[[1L]] <- quote(stats::model.frame) m <- eval.parent(m) Terms <- attr(m, "terms") X <- model.matrix(Terms, m, contrasts) cons <- attr(X, "contrasts") Xr <- qr(X)$rank Y <- model.response(m) if (!is.matrix(Y)) Y <- as.factor(Y) w <- model.weights(m) if (length(w) == 0L) if (is.matrix(Y)) w <- rep(1, dim(Y)[1L]) else w <- rep(1, length(Y)) lev <- levels(Y) if (is.factor(Y)) { counts <- table(Y) if (any(counts == 0L)) { empty <- lev[counts == 0L] warning(sprintf(ngettext(length(empty), "group %s is empty", "groups %s are empty"), paste(sQuote(empty), collapse = " ")), domain = NA) Y <- factor(Y, levels = lev[counts > 0L]) lev <- lev[counts > 0L] } if (length(lev) < 2L) stop("need two or more classes to fit a multinom model") if (length(lev) == 2L) Y <- as.integer(Y) - 1 else Y <- class.ind(Y) } if (summ == 1) { Z <- cbind(X, Y) z1 <- cumprod(apply(Z, 2L, max) + 1) Z1 <- apply(Z, 1L, function(x) sum(z1 * x)) oZ <- order(Z1) Z2 <- !duplicated(Z1[oZ]) oX <- (seq_along(Z1)[oZ])[Z2] X <- X[oX, , drop = FALSE] Y <- if (is.matrix(Y)) Y[oX, , drop = FALSE] else Y[oX] w <- diff(c(0, cumsum(w))[c(Z2, TRUE)]) print(dim(X)) } if (summ == 2) { Z <- summ2(cbind(X, Y), w) X <- Z$X[, 1L:ncol(X)] Y <- Z$X[, ncol(X) + 1L:ncol(Y), drop = FALSE] w <- Z$Y print(dim(X)) } if (summ == 3) { Z <- summ2(X, Y * w) X <- Z$X Y <- Z$Y[, 1L:ncol(Y), drop = FALSE] w <- rep(1, nrow(X)) print(dim(X)) } offset <- model.offset(m) r <- ncol(X) if (is.matrix(Y)) { p <- ncol(Y) sY <- Y %*% rep(1, p) if (any(sY == 0)) stop("some case has no observations") if (!censored) { Y <- Y/matrix(sY, nrow(Y), p) w <- w * sY } if (length(offset) > 1L) { if (ncol(offset) != p) stop("ncol(offset) is wrong") mask <- c(rep(FALSE, r + 1L + p), rep(c(FALSE, rep(TRUE, r), rep(FALSE, p)), p - 1L)) X <- cbind(X, offset) Wts <- as.vector(rbind(matrix(0, r + 1L, p), diag(p))) fit <- nnet(X, Y, w, Wts = Wts, mask = mask, size = 0, skip = TRUE, softmax = TRUE, censored = censored, rang = 0,MaxNWts=12000, ...) } else { mask <- c(rep(FALSE, r + 1L), rep(c(FALSE, rep(TRUE, r)), p - 1L)) fit <- nnet(X, Y, w, mask = mask, size = 0, skip = TRUE, softmax = TRUE, censored = censored, rang = 0,MaxNWts=12000, ...) } } else { if (length(offset) <= 1L) { mask <- c(FALSE, rep(TRUE, r)) fit <- nnet(X, Y, w, mask = mask, size = 0, skip = TRUE, entropy = TRUE, rang = 0,MaxNWts=12000, ...) } else { mask <- c(FALSE, rep(TRUE, r), FALSE) Wts <- c(rep(0, r + 1L), 1) X <- cbind(X, offset) fit <- nnet(X, Y, w, Wts = Wts, mask = mask, size = 0, skip = TRUE, entropy = TRUE, rang = 0, MaxNWts=12000, ...) } } fit$formula <- attr(Terms, "formula") fit$terms <- Terms fit$call <- call fit$weights <- w fit$lev <- lev fit$deviance <- 2 * fit$value fit$rank <- Xr edf <- ifelse(length(lev) == 2L, 1, length(lev) - 1) * Xr if (is.matrix(Y)) { edf <- (ncol(Y) - 1) * Xr if (length(dn <- colnames(Y)) > 0) fit$lab <- dn else fit$lab <- 1L:ncol(Y) } fit$coefnames <- colnames(X) fit$vcoefnames <- fit$coefnames[1L:r] fit$na.action <- attr(m, "na.action") fit$contrasts <- cons fit$xlevels <- .getXlevels(Terms, m) fit$edf <- edf fit$AIC <- fit$deviance + 2 * edf if (model) fit$model <- m class(fit) <- c("multinom", "nnet") if (Hess) fit$Hessian <- multinomHess(fit, X) fit}