FCM <- function(x, K, mybeta = 2, nstart = 1, iter_max = 100, eps = 1e-06) { ## FCM ## INPUTS ## x: input matrix n*d, n d-dim samples ## K: number of desired clusters ## Optional : ## mybeta : beta, exponent for u (defaut 2). ## nstart: how many random sets should be chosen(defaut 1) ## iter_max : The maximum number of iterations allowed. (default 100) ## ## OUTPUTS ## u: The fuzzy membership matrix = maxtrix of size n*K; ## g: matrix of size K*d of the centers of the clusters ## J: objective function ## histJ: all the objective function values in the iter process ## modified time: 2015-02-07 FCM_onetime <- function(x, init_centers, mybeta = 2, iter_max = 100, eps = 1e-06) { n = dim(x)[1] d = dim(x)[2] g = init_centers K = dim(g)[1] histJ = c() pasfini = 1 Jold = Inf D = matrix(0, n, K) for (j in 1:K) { D[, j] = rowSums(sweep(x, 2, g[j, ], "-")^2) } iter = 1 J_old = Inf while (pasfini) { s = (1/(D + eps))^(1/(mybeta - 1)) u = s/(s %*% matrix(1, K, K)) t1 = t(u^mybeta) %*% x t2 = t(u^mybeta) %*% matrix(1, n, d) V = t1/t2 g = V D = matrix(0, n, K) for (j in 1:K) { D[, j] = rowSums(sweep(x, 2, g[j, ], "-")^2) } J = sum(u^mybeta * D) pasfini = abs(J - Jold) > 0.001 && (iter < iter_max) Jold = J histJ = c(histJ, J) iter = iter + 1 } cluster_id = apply(u, 1, which.max) re = list(u, J, histJ, g, cluster_id) names(re) = c("u", "J", "histJ", "g", "cluster_id") return(re) } x = as.matrix(x) seeds = 1:nrow(x) id = sample(seeds, K) g = as.matrix(x[id, ]) re_best = FCM_onetime(x = x, init_centers = g, mybeta = mybeta, iter_max = iter_max, eps = eps) if (nstart > 1) { minJ = 0 i = 2 while (i <= nstart) { init_centers_id = sample(seeds, K) init_centers = as.matrix(x[init_centers_id, ]) run = FCM_onetime(x, init_centers = init_centers, mybeta = mybeta, iter_max = iter_max) if (run$J <= re_best$J) { re_best = run } i = i + 1 } } return(re_best) }
# 对于模糊聚类均值的公式及其推到，大致如下： #主要代码参见下面：（其中使用kmeans作比较。然后通过svm分类测验训练） # 设置伪随机种子 set.seed(100) # 生产数据样本 simple.data = function (n=200, nclass=2) { require(clusterGeneration) require(mvtnorm) # Center of Gaussians xpos = seq(-nclass*2, nclass*2, length=nclass) ypos = runif(nclass, min=-2*nclass, max=2*nclass) func = function(i,xpos,ypos,n) { # Create a random covariance matrix cov = genPositiveDefMat(2, covMethod="eigen", rangeVar=c(1, 10), lambdaLow=1, ratioLambda=10) # 保存随机数据 data = rmvnorm(n=n, mean=c(xpos[i], ypos[i]), sigma=cov$Sigma) # 保存每一次的结果 list(means=cbind(xpos[i], ypos[i]), covars=cov$Sigma, data=data,class=rep(i*1.0, n)) } # do call 合并列表 为 数据框 strL=do.call(rbind,lapply(1:nclass,func,xpos,ypos,n)) data=list() data$means=do.call(rbind,strL[,1]) data$covars = as.list(strL[,2]) data$data=do.call(rbind,strL[,3]) data$class=do.call(c,strL[,4]) # 返回 data } # 第一次随机产生u值 nr点个数 k 类别数 random.uij = function(k,nr) { # u = matrix(data=round(runif(k*nr,10,20)),nrow=k,ncol=nr, dimnames=list(paste(\'u\',1:k,sep=""),paste(\'x\',1:nr,sep=\'\'))) tempu = function(x) { ret = round(x/sum(x),4) # 保证每一列之和为1 ret[1] = 1-sum(ret[-1]) ret } apply(u,2,tempu) } # 计算 点矩阵 到 中心的距离 dist_cc_dd = function(cc,dd) { # cc 为 中心点 dd 为样本点值 temp = function(cc,dd) { # 计算每一个中心点与每一个点的距离 temp1 = function(index) { sqrt(sum(index^2)) } # 结果向量以列存放，后面的结果需要转置，按行存储 apply(dd-cc,2,temp1) } # 将结果转置 t(apply(cc,1,temp,dd)) } # 模糊均值聚类 fuzzy.cmeans = function(data,u,m=3) { # 简单的判断，可以不要 if (is.array(data) || is.matrix(data)) { data = as.data.frame(data) } # nr = nrow(data) # nc = ncol(data) # while (J > lim && step < steps) # { # step = step + 1 # uij 的 m 次幂 um = u^m rowsum = apply(um,1,sum) # 求中心点 ci cc = as.matrix(um/rowsum) %*% as.matrix(data) # rownames(cc)=paste(\'c\',1:k,sep=\'\') # colnames(cc)=paste(\'x\',1:nc,sep=\'\') # 计算 J 值 distance = dist_cc_dd(cc,t(data)) J = sum(distance^2 * um) # cc_temp = matrix(rep(cc,each=nr),ncol=2) # dd_temp = NULL # lapply(1:k,function(i){dd_temp <<- rbind(dd_temp,data)}) # dist = apply((dd_temp-cc_temp)^2,1,sum) # um_temp = as.vector(t(um)) # J = um_temp %*% dist # 计算幂次系数，后面需要使用m != 1 t = -2 / (m - 1) # 根据公式 计算 tmp = distance^t colsum = apply(tmp,2,sum) mat = rep(1,nrow(cc)) %*% t(colsum) # 计算 uij，如此u的每一列之和为0 u = tmp / mat # } # u # 保存一次迭代的结果值 list(U = u,C = cc,J = J) } # 设置初始化参数 n = 100 k = 4 dat = simple.data(n,k) nr = nrow(dat$data) m = 3 limit = 1e-4 max_itr=50 # 随机初始化 uij u = random.uij(k,nr) results = list() data=dat$data # 迭代计算 收敛值 for (i in 1 : max_itr) { results[[i]] = fuzzy.cmeans(dat$data,u,m) if (i != 1 && abs((results[[i]]$J - results[[i-1]]$J)) < limit) { break } u = results[[i]]$U } # 做散点图 require(ggplot2) data=as.data.frame(dat$data,stringsAsFactors=FALSE) data=cbind(data,dat$class) nc = ncol(data) colnames(data)=paste(\'x\',1:nc,sep=\'\') # par(mar=rep(2,4)) p=ggplot(data,aes(x=x1,y=x2,color=factor(x3))) p+geom_point()+xlab(\'x轴\')+ylab(\'y轴\')+ggtitle(\'scatter points\') # plot(dat$data,col=factor(dat$class)) # points(results[[i]]$C,pch=19,col=1:uniqe(dat$class)) # Sys.sleep(1) # 计算聚类与原始类的误差比率 tclass=apply(results[[i]]$U,2,function(x){which(x==max(x))}) tclass[tclass==2]=5 tclass[tclass==3]=6 tclass[tclass==4]=7 tclass[tclass==5]=4 tclass[tclass==6]=2 tclass[tclass==7]=3 freq=table(dat$class,tclass) (sum(freq)-sum(diag(freq))) / sum(freq) # 训练 svm model svm_test = function() { library(e1071) svm.fit = svm(dat$data,dat$class) r.fit = predict(svm.fit, dat$data) diff.class = round(as.numeric(r.fit)) - as.numeric(dat$class) i.misclass = which(abs(diff.class) > 0) n.misclass = length(i.misclass) f.misclass = n.misclass/length(dat$class) } # 同一数据，使用 kmeans 聚类 kmeans_test = function() { k.fit = kmeans(x=dat$data,4) tclass=k.fit$cluster tclass[tclass==2]=5 tclass[tclass==3]=6 tclass[tclass==4]=7 tclass[tclass==5]=3 tclass[tclass==6]=4 tclass[tclass==7]=2 freq=table(dat$class,tclass) (sum(freq)-sum(diag(freq))) / sum(freq) } # kmeans 和 fuzzy c means
请发表评论