40 付録2：混合正規分布に適用するEMアルゴリズムのRソースコード #read data from file in *.csv from the path string options(warn=-1) data <- read.csv("D://result//case2//sample.csv",header=T) #working days of TMC data <- subset(data, working == 1 & peak == 1) #data <- subset(data, working == 1 & peak == 0) #No working days of TMC #data <- subset(data, working == 0 & peak == 1) #data <- subset(data, working == 0 & peak == 0) a <- c(0.3,0.7) u <- c(397.97,360.71) sigma2 <- c(199.40*199.40,131.73*131.73) GaussCluster <- function(data,k,a,u,sigma2) { #data:save bus travel time which will be separated #k:number of Gaussian Mixed Distribution #a0: priorior distribution and related vectors #u0:initial values of means #sigma2:variance of travel time #https://blog.csdn.net/qq_29737811/article/details/78836427 mat = matrix(rep(data,k),ncol = k) #Length of sample data N = length(data) i=0 while(TRUE) { u0 <- u # sigma <- sqrt(sigma2) sigma0 <- sigma umatrix <- matrix(rep(u,N),nrow=N,byrow = T) p <- t(sapply(data,dnorm,u,sigma)) amatrix <- matrix(rep(a,N),nrow=N,byrow = T) r <- (p * amatrix) / rowSums(p * amatrix) u <- colSums(data*r)/colSums(r) sigma2 <- colSums(r*(mat-umatrix)^2)/colSums(r) a <- colSums(r)/length(data) sumU <- sum((u-u0)^2) sumsigma <- sum((sigma-sigma0)^2) if((sumU <= 1e-4) & (sumsigma <= 1e-4)) { break } } cluster <- which(r == apply(r,1,max),arr.ind = T) return(list(u = u,d = sigma2,a = a,cluster = cluster)) } results <- GaussCluster(data[,6],2,a,u,sigma2) cat("Mean of Travel Time:",round(results$u, digits = 3),"¥n") cat("Standard Deviation of Travel Time:",round(sqrt(results$d), digits = 3),"¥n") cat("Composite Ratio:",round(results$a, digits = 3),"¥n")
元のページ ../index.html#46