SamNum = 100; % 總樣本數
TestSamNum = 101; % 測試樣本數
InDim = 1; % 樣本輸入維數
ClusterNum = 10; % 隱節點數，即聚類樣本數
Overlap = 1.0; % 隱節點重疊系數

% 根據目標函數獲得樣本輸入輸出
rand('state',sum(100*clock))
NoiseVar = 0.1;
Noise = NoiseVar*randn(1,SamNum);
SamIn = 8*rand(1,SamNum)-4;
SamOutNoNoise = 1.1*(1-SamIn+2*SamIn.^2).*exp(-SamIn.^2/2);
SamOut = SamOutNoNoise + Noise;

TestSamIn = -4:0.08:4;
TestSamOut = 1.1*(1-TestSamIn+2*TestSamIn.^2).*exp(-TestSamIn.^2/2);

figure;
hold on;
grid;
plot(SamIn,SamOut,'k+')
plot(TestSamIn,TestSamOut,'k--')
xlabel('Input x');
ylabel('Output y');

Centers = SamIn(:,1:ClusterNum);

NumberInClusters = zeros(ClusterNum,1); % 各類中的樣本數，初始化為零
IndexInClusters = zeros(ClusterNum,SamNum); % 各類所含樣本的索引號
flag=1;
while flag
NumberInClusters = zeros(ClusterNum,1); % 各類中的樣本數，初始化為零
IndexInClusters = zeros(ClusterNum,SamNum); % 各類所含樣本的索引號

% 按最小距離原則對所有樣本進行分類
for i = 1:SamNum
    AllDistance = dist(Centers',SamIn(:,i));
    [MinDist,Pos] = min(AllDistance);
    NumberInClusters(Pos) = NumberInClusters(Pos) + 1;
    IndexInClusters(Pos,NumberInClusters(Pos)) = i;
end
% 保存舊的聚類中心
OldCenters = Centers;

for i = 1:ClusterNum
    Index = IndexInClusters(i,1:NumberInClusters(i));
    Centers(:,i) = mean(SamIn(:,Index)')';
end
% 判斷新舊聚類中心是否一致，是則結束聚類
EqualNum = sum(sum(Centers==OldCenters));
if EqualNum == InDim*ClusterNum
break
end
end

% 計算各隱節點的擴展常數（寬度）
AllDistances = dist(Centers',Centers); % 計算隱節點數據中心間的距離（矩陣）
Maximum = max(max(AllDistances)); % 找出其中最大的一個距離
for i = 1:ClusterNum % 將對角線上的0 替換為較大的值
AllDistances(i,i) = Maximum+1;
end
Spreads = Overlap*min(AllDistances)'; % 以隱節點間的最小距離作為擴展常數

% 計算各隱節點的輸出權值
Distance = dist(Centers',SamIn); % 計算各樣本輸入離各數據中心的距離
SpreadsMat = repmat(Spreads,1,SamNum);
HiddenUnitOut = radbas(Distance./SpreadsMat); % 計算隱節點輸出陣
HiddenUnitOutEx = [HiddenUnitOut' ones(SamNum,1)]'; % 考慮偏移
W2Ex = SamOut*pinv(HiddenUnitOutEx); % 求廣義輸出權值
W2 = W2Ex(:,1:ClusterNum); % 輸出權值
B2 = W2Ex(:,ClusterNum+1); % 偏移

% 測試
TestDistance = dist(Centers',TestSamIn);
TestSpreadsMat = repmat(Spreads,1,TestSamNum);
TestHiddenUnitOut = radbas(TestDistance./TestSpreadsMat);
TestNNOut = W2*TestHiddenUnitOut+B2;
plot(TestSamIn,TestNNOut,'k-');
W2;
B2;