1. 基本的なFCM実装

function [U, centers, J_hist] = fuzzyCMeans(data, k, fuzz, maxIter, tol)
    % data: n×d 行列
    % k: クラスタ数
    % fuzz: ファジィ度（通常1.5～2.5）
    % maxIter: 最大反復回数
    % tol: 収束閾値
    
    [n, d] = size(data);
    
    % ランダム初期化
    U = rand(n, k);
    U = U ./ sum(U, 2);
    
    J_hist = zeros(maxIter, 1);
    
    for t = 1:maxIter
        % クラスタ中心更新
        um = U .^ fuzz;
        centers = (um' * data) ./ sum(um)' * ones(1, d);
        
        % 距離行列
        D = zeros(n, k);
        for i = 1:n
            for j = 1:k
                D(i,j) = norm(data(i,:) - centers(j,:));
            end
        end
        
        % 所属度更新
        D = max(D, 1e-6); % ゼロ除算回避
        U_prev = U;
        U = 1 ./ (D .^ (2/(fuzz-1)) .* sum(1./(D .^ (2/(fuzz-1))), 2));
        
        % 収束判定
        if norm(U - U_prev, 'fro') < tol
            J_hist = J_hist(1:t);
            break;
        end
        
        % 目的関数値保存
        J_hist(t) = sum(sum((U.^fuzz) .* (D.^2)));
    end
end

2. カーネルFCM（高斯カーネル版）

function [U, centers, J_hist] = kernelFCM(data, k, fuzz, maxIter, tol, sigma)
    [n, ~] = size(data);
    
    % カーネル行列
    K = exp(-squareform(pdist(data)).^2 / (2*sigma^2));
    
    % 初期化
    U = rand(n, k);
    U = U ./ sum(U, 2);
    
    J_hist = zeros(maxIter, 1);
    
    for t = 1:maxIter
        % カーネル空間での中心
        um = U .^ fuzz;
        centers = (um' * K) ./ sum(um)';
        
        % カーネル距離
        distK = zeros(n, k);
        for i = 1:n
            for j = 1:k
                distK(i,j) = K(i,i) - 2*centers(j,:)*K(:,i) + centers(j,:)*K*centers(j,:)';
            end
        end
        
        % 所属度更新
        U_prev = U;
        U = 1 ./ (distK .^ (1/(fuzz-1)) .* sum(1./(distK .^ (1/(fuzz-1))), 2));
        
        if norm(U - U_prev, 'fro') < tol
            J_hist = J_hist(1:t);
            break;
        end
        
        J_hist(t) = sum(sum((U.^fuzz) .* distK));
    end
end

3. 実験用データ生成と実行例

% 2次元ガウス分布混合
rng(2024);
mu1 = [2, 3];  mu2 = [7, 8];
cov1 = [1 0.5; 0.5 1];  cov2 = [1 -0.3; -0.3 1];
X = [mvnrnd(mu1, cov1, 150); mvnrnd(mu2, cov2, 150)];

% パラメータ
k = 2; fuzz = 2; maxIter = 150; tol = 1e-4; sigma = 1.5;

% FCM実行
tic; [U_fcm, C_fcm, J_fcm] = fuzzyCMeans(X, k, fuzz, maxIter, toc); t_fcm = toc;

% KFCM実行
tic; [U_kfcm, C_kfcm, J_kfcm] = kernelFCM(X, k, fuzz, maxIter, tol, sigma); t_kfcm = toc;

4. 結果可視化

figure('Position', [100 100 900 400]);

subplot(1,2,1);
scatter(X(:,1), X(:,2), 30, max(U_fcm,[],2), 'filled');
hold on;
plot(C_fcm(:,1), C_fcm(:,2), 'kp', 'MarkerSize', 15, 'LineWidth', 2);
title(sprintf('FCM\n反復: %d, 時間: %.3fs', numel(J_fcm), t_fcm));
colormap(gca, lines(k));

subplot(1,2,2);
scatter(X(:,1), X(:,2), 30, max(U_kfcm,[],2), 'filled');
hold on;
plot(C_kfcm(:,1), C_kfcm(:,2), 'kp', 'MarkerSize', 15, 'LineWidth', 2);
title(sprintf('KFCM (σ=%.1f)\n反復: %d, 時間: %.3fs', sigma, numel(J_kfcm), t_kfcm));
colormap(gca, lines(k));

5. 高速化と応用例

5.1 K-means初期化による高速化

function U = initWithKmeans(data, k)
    [~, idx] = kmeans(data, k, 'Replicates', 10);
    U = full(sparse(1:numel(idx), idx, 1));
end

5.2 画像セグメンテーション

% グレー画像読み込み
img = imread('rice.png');
pixels = double(img(:));

% KFCMセグメンテーション
[U, ~, ~] = kernelFCM(pixels, 3, 2, 100, 1e-4, 30);

% セグメント画像生成
seg = reshape(U, size(img));
segmented = ind2rgb(uint8(255*max(seg,[],3)), jet(3));
imshow(segmented);

5.3 時系列異常検知

% 振動信号
load('bearing.mat');
features = [mean(sig), std(sig), kurtosis(sig), skewness(sig)];

% 2クラスタFCM
[U, ~, ~] = fuzzyCMeans(features, 2, 2, 100, 1e-4);
anomaly = find(U(:,2) > 0.8);

6. パフォーマンス比較表