R語言 模糊c均值(FCM)算法程序（轉）

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