function [ c, sigma , W_output ] = SOFNN( X, d, Kd )
%SOFNN Self-Organizing Fuzzy Neural Networks
%Input Parameters
% X(r,n) - rth traning data from nth observation
% d(n) - the desired output of the network (must be a row vector)
% Kd(r) - predefined distance threshold for the rth input

%Output Parameters
% c(IndexInputVariable,IndexNeuron)
% sigma(IndexInputVariable,IndexNeuron)
% W_output is a vector

%Setting up Parameters for SOFNN
SigmaZero=4;
delta=0.12;
threshold=0.1354;
k_sigma=1.12;

%For more accurate results uncomment the following
%format long;

%Implementation of a SOFNN model
[size_R,size_N]=size(X);
%size_R - the number of input variables
c=[];
sigma=[];
W_output=[];
u=0; % the number of neurons in the structure
Q=[];
O=[];
Psi=[];
for n=1:size_N
    x=X(:,n);   
    if u==0 % No neuron in the structure?
        c=x;
        sigma=SigmaZero*ones(size_R,1);
        u=1;
        Psi=GetMePsi(X,c,sigma);
        [Q,O] = UpdateStructure(X,Psi,d);
        pT_n=GetMeGreatPsi(x,Psi(n,:))';        
    else
        [Q,O,pT_n] = UpdateStructureRecursively(X,Psi,Q,O,d,n);
    end;

    KeepSpinning=true;
    while KeepSpinning
        %Calculate the error and if-part criteria
        ae=abs(d(n)-pT_n*O); %approximation error
        [phi,~]=GetMePhi(x,c,sigma);
        [maxphi,maxindex]=max(phi); % maxindex refers to the neuron's index

        if ae>delta
            if maxphi<threshold
                %enlarge width
                [minsigma,minindex]=min(sigma(:,maxindex));
                sigma(minindex,maxindex)=k_sigma*minsigma;
                Psi=GetMePsi(X,c,sigma);
                [Q,O] = UpdateStructure(X,Psi,d);
                pT_n=GetMeGreatPsi(x,Psi(n,:))';                              
            else
                %Add a new neuron and update structure
                ctemp=[];
                sigmatemp=[];
                dist=0;
                for r=1:size_R
                    dist=abs(x(r)-c(r,1));
                    distIndex=1;
                    for j=2:u
                        if abs(x(r)-c(r,j))<dist
                            distIndex=j;
                            dist=abs(x(r)-c(r,j));
                        end;
                    end;
                    if dist<=Kd(r)
                        ctemp=[ctemp; c(r,distIndex)];
                        sigmatemp=[sigmatemp ; sigma(r,distIndex)];
                    else
                        ctemp=[ctemp; x(r)];
                        sigmatemp=[sigmatemp ; dist];
                    end;
                end;
                c=[c ctemp];
                sigma=[sigma sigmatemp];
                Psi=GetMePsi(X,c,sigma);
                [Q,O] = UpdateStructure(X,Psi,d);
                KeepSpinning=false;
                u=u+1;
            end;
        else
            if maxphi<threshold
                %enlarge width
                [minsigma,minindex]=min(sigma(:,maxindex));
                sigma(minindex,maxindex)=k_sigma*minsigma;
                Psi=GetMePsi(X,c,sigma);
                [Q,O] = UpdateStructure(X,Psi,d);
                pT_n=GetMeGreatPsi(x,Psi(n,:))';               
            else
                %Do nothing and exit the while
                KeepSpinning=false;               
            end;
        end;        
    end;
end;
W_output=O;
end

function [Q_next, O_next,pT_n] = UpdateStructureRecursively(X,Psi,Q,O,d,n)
%O=O(t-1) O_next=O(t)
p_n=GetMeGreatPsi(X(:,n),Psi(n,:));
pT_n=p_n';
ee=abs(d(n)-pT_n*O); %|e(t)|
temp=1+pT_n*Q*p_n;
ae=abs(ee/temp);

if ee>=ae
    L=Q*p_n*(temp)^(-1);
    Q_next=(eye(length(Q))-L*pT_n)*Q;
    O_next=O + L*ee;
else
    Q_next=eye(length(Q))*Q;
    O_next=O;
end;
end

function [ Q , O ] = UpdateStructure(X,Psi,d)
GreatPsiBig = GetMeGreatPsi(X,Psi);

%M=u*(r+1)
%n - the number of observations
[M,~]=size(GreatPsiBig);

%Others Ways of getting Q=[P^T(t)*P(t)]^-1
%**************************************************************************
%opts.SYM = true;
%Q = linsolve(GreatPsiBig*GreatPsiBig',eye(M),opts);
%
%Q = inv(GreatPsiBig*GreatPsiBig');
%Q = pinv(GreatPsiBig*GreatPsiBig');
%**************************************************************************
Y=GreatPsiBig\eye(M);
Q=GreatPsiBig'\Y;

O=Q*GreatPsiBig*d';
end


%This function works too with x
% (X=X and Psi is a Matrix) - Gets you the whole GreatPsi
% (X=x and Psi is the row related to x) - Gets you just the column related with the observation
function [GreatPsi] = GetMeGreatPsi(X,Psi)
%Psi - In a row you go through the neurons and in a column you go through number of
%observations **** Psi(#obs,IndexNeuron) ****
GreatPsi=[];
[N,U]=size(Psi);
for n=1:N
    x=X(:,n);
    GreatPsiCol=[];
    for u=1:U
        GreatPsiCol=[ GreatPsiCol ; Psi(n,u)*[1; x] ];
    end;
    GreatPsi=[GreatPsi GreatPsiCol];
end;
end


function [phi, SumPhi]=GetMePhi(x,c,sigma)
[r,u]=size(c);
%u - the number of neurons in the structure
%r - the number of input variables

phi=[];
SumPhi=0;

for j=1:u % moving through the neurons
    S=0;
    for i=1:r % moving through the input variables
        S = S + ((x(i) - c(i,j))^2) / (2*sigma(i,j)^2);
    end;   
    phi = [phi exp(-S)];        
    SumPhi = SumPhi + phi(j);   %phi(u)=exp(-S)
end;
end


%This function works too with x, it will give you the row related to x
function [Psi] = GetMePsi(X,c,sigma)
[~,u]=size(c);
[~,size_N]=size(X);
%u - the number of neurons in the structure
%size_N - the number of observations
Psi=[];
for n=1:size_N        
    [phi, SumPhi]=GetMePhi(X(:,n),c,sigma);   
    PsiTemp=[];
    for j=1:u
        %PsiTemp is a row vector ex: [1 2 3]
        PsiTemp(j)=phi(j)/SumPhi;
    end;

    Psi=[Psi; PsiTemp];   
    %Psi - In a row you go through the neurons and in a column you go through number of
    %observations **** Psi(#obs,IndexNeuron) ****
end;

end

自己调试一下

嗯  谢谢你热心的帮助，我初学matlab和神经网络这块，这个代码可以执行，但怎样进行训练和预测可以简单给我说说吗？

嗯  谢谢你热心的帮助，我初学matlab和神经网络这块，这个代码可以执行，但怎样进行训练和预测可以简单给我说说吗？