シンプレックス法を行ったlp.assign関数とlpSolveパッケージ。
lp.assign function (cost.mat, direction = "min", presolve = 0, compute.sens = 0) { if (!is.matrix(cost.mat)) stop("Matrix of costs required.") if (is.data.frame(cost.mat)) cost.mat <- as.matrix(cost.mat) nr <- nrow(cost.mat) nc <- ncol(cost.mat) rnum.signs <- rep(3, nr) row.rhs <- rep(1, nr) cnum.signs <- rep(3, nc) col.rhs <- rep(1, nc) if (direction == "min") direction <- as.integer(0) else if (direction == "max") direction <- as.integer(1) else stop("Direction must be 'min' or 'max'") varcount <- as.integer(nr * nc) objective <- as.double(c(0, c(t(cost.mat)))) const.count <- as.integer(nr + nc) intcount <- as.integer(varcount) intvec <- as.integer(1:varcount) objval <- as.double(0) int.count <- nc * nr integers <- as.integer(numeric(int.count)) solution <- as.double(numeric(nc * nr)) status <- as.integer(0) sens.coef.from <- sens.coef.to <- 0 duals <- duals.from <- duals.to <- 0 if (compute.sens) { sens.coef.from <- sens.coef.to <- numeric(varcount) duals <- duals.from <- duals.to <- numeric(varcount + const.count) } lps.out <- .C("lp_transbig", direction = direction, rcount = as.integer(nr), ccount = as.integer(nc), costs = objective, rsigns = as.integer(rnum.signs), rrhs = as.double(row.rhs), csigns = as.integer(cnum.signs), crhs = as.double(col.rhs), objval = objval, int.count = int.count, integers = integers, solution = solution, presolve = as.integer(presolve), compute.sens = as.integer(compute.sens), sens.coef.from = as.double(sens.coef.from), sens.coef.to = as.double(sens.coef.to), duals = as.double(duals), duals.from = as.double(duals.from), duals.to = as.double(duals.to), status = status, PACKAGE = "lpSolve") lps.out$solution = matrix(lps.out$solution, nr, nc, byrow = TRUE) if (length(duals) > 0) { lps.out$sens.coef.from <- matrix(lps.out$sens.coef.from, nr, nc, byrow = TRUE) lps.out$sens.coef.to <- matrix(lps.out$sens.coef.to, nr, nc, byrow = TRUE) lps.out$duals <- matrix(lps.out$duals, nr, nc, byrow = TRUE) lps.out$duals.from <- matrix(lps.out$duals.from, nr, nc, byrow = TRUE) lps.out$duals.to <- matrix(lps.out$duals.to, nr, nc, byrow = TRUE) } lps.out$costs <- cost.mat if (any(names(version) == "language")) class(lps.out) <- "lp" else oldClass(lps.out) <- "lp" lps.out } <environment: namespace:lpSolve>
データの準備
studentid <- c("RA", "NA", "MA", "MA", "YA", "HA", "TI", "YI", "KI", "HI", "YI", "YI", "YI", "EI", "KI", "YI", "HI", "YU", "YO", "KO", "KO", "YO", "SO", "YO", "JO", "TO", "NO", "SO", "YO", "TO", "MO", "TK", "TK", "RK", "KK", "YK", "RK", "MK", "YK", "KG", "KS", "YS", "KS", "KS", "YS", "YS", "ES", "KS", "SS", "MD", "KT", "AT", "YT", "ST", "TN", "TN", "KN", "YN", "TN", "MN", "RN", "KH", "YH", "AH", "YH", "MH", "GF", "YF", "HF", "MF", "RF", "KH", "TM", "TM", "KM", "NM", "SM", "HM", "AM", "YM", "HM", "YM", "NM", "TY", "TY", "NY", "AY", "TY", "FY", "NY", "RY", "AY", "KW", "KY", "TA", "MO", "HO", "YO", "HK", "HK", "NS", "MS", "KT", "TH", "JF", "HH", "YK", "KD", "HF", "TY") studentid <- paste(1:length(studentid), studentid, sep="") finalexpect <- as.matrix(read.csv("simplexexpectdata.csv")) sscourse <- as.matrix(read.csv("sscourse.csv"))[, -1] ssgroup <- as.matrix(read.csv("ssgroup.csv"))[, -1] ssexpect <- as.matrix(read.csv("ssexpect.csv"))[, -1] niter <- ncol(sscourse) #上位者 highscorer <- c(74, 27, 102, 92, 42, 7, 80, 72, 35, 94, 81, 46) #行が学生、列がシミュレーション回数、割当コース、希望、小グループ simplexarray <- array(0, c(110, niter, 3)) dimnames(simplexarray) <- list(studentid, 1:niter, c("course", "expect", "group")) simplexarray[ highscorer, , "course"] <- c(10, 14, 16, 21, 24, 38, 47, 75, 77, 79, 98, 101) simplexarray[-highscorer, , "course"] <- sscourse simplexarray[ highscorer, , "expect"] <- 3 simplexarray[-highscorer, , "expect"] <- ssexpect simplexarray[ highscorer, , "group"] <- c(2, 3, 3, 4, 4, 7, 9, 14, 14, 15, 18, 19) simplexarray[-highscorer, , "group"] <- ssgroup
各コースの希望の累計数
count.expect <- function(vec){ exA <- sum(vec == 3) exB <- sum(vec == 2) exD <- sum(vec == -15) return(c(exA, exB, exD)) } expect_res <- apply(finalexpect, 2, count.expect) #png("No_course_expect.png", 1200, 300) barplot(expect_res, col=c("red", "yellow", "green"), main="No of ") xaxispoint <- barplot(expect_res, plot=FALSE) axis(1, at=xaxispoint, labels=1:ncol(expect_res), tick=FALSE, pos=0.03, las=2, cex.axis=0.8) #dev.off()
コースが確定してしまう人
#コースが確定してしまう人 d <- simplexarray[which(sapply(apply(simplexarray[, , "course"], 1, unique), length) == 1), 1, "course"]
どれくらいのシミュレーション回数で収束するか確認。
もっとも取りうるコース数が多い108番で確認。
#asymp.array <- array(0, c(niter, dim(simplexarray)[1], length(studentid))) #これはメモリ的に無理 asymp.array <- array(0, c(niter, dim(simplexarray)[1], 1)) for(who in 1){ tempmat <- matrix(0, niter, dim(simplexarray)[1]) for(ni in 1:niter){ temp <- simplexarray[c(108)[who], ,"course"] tempmat[ni, temp[ni]] <- 1 print(ni) } asymp.array[, , who] <- apply(tempmat, 2, cumsum) / (1:niter) matplot(asymp.array[, , who], type="l", xlab="No. simulation trials", ylab="probability") title(paste("course assign probability simulation")) }
どのコースになるか、確率の収束
prob.course <- diag(0, dim(simplexarray)[1]) dimnames(prob.course) <- list(studentid, 1:dim(simplexarray)[1]) for(st in 1:dim(simplexarray)[1]){ coursecountvec <- numeric(dim(simplexarray)[1]) for(cs in 1:1:dim(simplexarray)[1]){ coursecountvec[cs] <- sum(simplexarray[st, , "course"] == cs) / niter print(c(st, cs)) } prob.course[st, ] <- coursecountvec } #write.csv(prob.course, file="prob.course.csv")
各個人が属するグループの確率
indvgroup <- matrix(0, nrow(simplexarray[, , "group"]), length(unique(c(simplexarray[, , "group"])))) dimnames(indvgroup) <- list(studentid, 1:length(unique(c(simplexarray[, , "group"])))) niter <- dim(simplexarray)[2] for(st in 1:nrow(simplexarray[, , "group"])){ for(gr in 1:length(unique(c(simplexarray[, , "group"])))){ indvgroup[st, gr] <- sum(simplexarray[st, , "group"] == gr) / niter #print(c(st, gr)) } }
二人選んで、一緒のグループになる確率
partnerres <- diag(0, nrow(simplexarray[, , "group"])) dimnames(partnerres) <- list(studentid, studentid) niter <- dim(simplexarray)[2] for(g in 1:nrow(simplexarray[, , "group"])){ partnerres[, g] <- t(apply(t(t(simplexarray[, , "group"]) == simplexarray[g, , "group"]), 1, cumsum) / 1:niter)[, ncol(simplexarray[, , "group"])] #print(g) }
A希望になる確率
expect.count <- function(vec, expect.score){ return(sum(vec == expect.score)) } resA <- apply(simplexarray[, , "expect"], 1, expect.count, expect.score=3) resA <- resA / niter #png("expect_probability.png", 1200, 300) barplot(resA, col=rainbow(length(resA)), axisname=FALSE) xaxispoint <- barplot(resA, plot=FALSE) axis(1, at=xaxispoint, labels=rownames(simplexarray), tick=FALSE, pos=0.01, las=2, cex.axis=0.8) #dev.off()
A希望の数と割当結果
NoA <- apply(finalexpect, 1, count.expect)[1,] Acoef <- cbind(NoA[-highscorer], resA[-highscorer]) Alist <- lapply(mapply(rep, 0, 4:6), unique) for(i in 4:6){ if(i != 6){ Alist[[i - 3]] <- Acoef[which(Acoef[, 1] == i), 2] }else{ Alist[[i - 3]] <- Acoef[which(Acoef[, 1] >= i), 2] } } Ag <- rep(1:length(sapply(Alist, length)), sapply(Alist, length)) pairwise.t.test(unlist(Alist), Ag, p.adj="bonf")
誰と同じになるかヒートマップ
library(gplots) #png("heatmap_twopair.png", 1200, 1200) heatmap.2(partnerres, scale="none", col=c(colorpanel(100, low="white", high="black"), "red"), symkey=FALSE, key=TRUE, keysize=0.3, density.info="none", tracecol=NULL, cexRow=1, cexCol=1, distfun=function(a){dist(a, method="maximum")}, hclustfun=function(b){hclust(b, method="complete")}) #dev.off() #png("heatmap_group.png", 300, 1200) heatmap.2(indvgroup, scale="none", col=c(colorpanel(100, low="white", high="black"), "red"), symkey=FALSE, key=FALSE, keysize=0.2, density.info="none", trace="none", dendrogram="none", cexRow=1, cexCol=1.5, distfun=function(a){dist(a, method="maximum")}, hclustfun=function(b){hclust(b, method="complete")}) #dev.off()
