%%
clear
clc
tic,
%[x,y]=data;
x=[1 1 1;
1 2 3];
y=[2 3 4]; %%%%%--數據顯示，輸入為－兩輸入，輸出為－單輸出。--------樣本為p2組
[p1,p2]=size(x);
%－ 一。首先要對樣本進行聚類分析，以此來確定模糊規則個數。利用K-means法對樣本聚類。

%????此處的K－ means 法待加

%－ 二。建立模糊推理系統

% 隸屬度函數個數-－模糊規則個數
k=5;
% 初始化四個隸屬度函數的參數A,B及輸出層初始權值W

for i=1:p1; 
for j=1:k;
m(i,j)=1+0.6*rands(1);
b(i,j)=1+0.6*rands(1);
end
end
for j=1:k;
w(j)=100+rand(1);
end
%%%---推理計算輸出值
for q=1:p2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%-----用同一隸屬度參數對 輸入樣本 X 累計計算
% 選用高斯函數作為隸屬度，求隸屬度，共 size(x,2)+k 個。x(1) K個，x(2) K個
for i=1:p1;
for j=1:k;
u(i,j)=gaussmf(x(i,q),[m(i,j),b(i,j)]);
end
end
% 模糊推理計算：a21,a22.幾個隸屬度函數，得出幾個值，本例中兩個.
%%%%----由以前的取小做法改為相乘—prod(x,1) or prod(x,2)———
for i=1:k;
v(i)=1;
j=1;
while j<=p1;
v(i)=v(i)*u(j,i);
j=j+1;
end
end
% 歸一化計算模糊推理的值；相當于已經除去了經典去模糊輸出的分母值
%for i=1:k;
%a3(i)=a2(i)/sum(a2);
%end
% 系統輸出
out1(q)=w*v';
e(q)=(y(q)-out1(q));
end
out=out1

%－ 三。參數修正過程。 增加方式，非批處理方式迭代
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%－－－－－－－－---------------------誤差反向傳播過程－－－-----------------------------------------
% 取誤差函數：E＝(1/2)*sumsqr(t-y)
E=(1/2)*sumsqr(y-out)
EE=E;
e
% e=sum(y-out)
lr=0.3; % c2=zeros(2,2);
%%%%----------------------------------------誤差反傳后的參數修正過程-------------------
r=1;
p=1;
s=1000;
while p<=s & EE>1e-10
%%%%%%%%%%%%%_____隸屬度參數 M. B 輸出層權值參數 W 的修正過程_____%%%%%%%%%%%%
%%1.--W
wc=zeros(1,k);
for i=1:k;
wc(i)=lr*e(r)*v(i);
end
%%2.--M
mc=zeros(p1,k);
for i=1:p1;
for j=1:k;
mc(i,j)=2*lr*e(r) * w(j) * (v(j)/u(i,j)) * exp(-((x(i,r)-m(i,j)).^2)/(b(i,j).^2))* (x(i,r)-m(i,j))/(b(i,j).^2);
end
end
%%3.--B
bc=zeros(p1,k);
for i=1:p1;
for j=1:k;
bc(i,j)=2*lr*e(r)* w(j) * (v(j)/u(i,j)) * exp(-((x(i,r)-m(i,j)).^2)/(b(i,j).^2)) * ((x(i,r)-m(i,j)).^2)/(b(i,j).^3);
end
end
% 4.參數修正 m b w
m=m+mc;
b=b+bc;
w=w+wc;
%%%%%%%%%%%_______利用修正后的參數重新計算_____________%%%%%%%%%%%%%%%%%%%%%
% 5.利用修正過的參數重新計算輸出
for q=1:p2; 
for i=1:p1;
for j=1:k;
u(i,j)=gaussmf(x(i,q),[m(i,j),b(i,j)]);
end
end
for i=1:k;
v(i)=1;
j=1;
while j<=p1;
v(i)=v(i)*u(j,i);
j=j+1;
end
end
out1(q)=w*v';
end
out=out1;
p=p+1;
EE=(1/2)*sumsqr(y-out);
E(p)=EE;
r=r+1;
if r>p2
r=1;
end
e(r)=(y(r)-out(r));
end
%%%%%%%%%%%%%%%%%%%________________當誤差或迭代步數滿足要求后得到結果_________________%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m,b,w,E_out=EE,e
epoch=1:size(E,2);
figure
plot(epoch,E,'-r');
axis([0 1.5*s min(E) max(E)]);
set(gca,'fontsize',8);
set(gca,'xtick',0:s/10:1.5*s);
%set(gca,'ytick',1e-30:1e5:1e5);
%set(gcf,'color','b')
title('誤差變化曲線');xlabel('步數');ylabel('誤差');
toc