function output = DRNNX(Is,Os)
  % Is 是特征矩陣，大小為 n * 8， n為樣本個數, 8 為特征向量的長度
  % Os 為輸出參考矩陣, 大小為 n* 1
  % 輸出為cell變量{Woij,Wdij,WIij},Woij = output{1},Wdij = output{2}, WIij = output{3}

NumI = 8;                                                                  % i = 8;
NumJ = 17;                                                                 % j = 17;
Wi = ones(NumI,NumJ);                                                      % Wij;
X = zeros(size(Is,1),NumJ);                                                % Xij;
Wd = ones(1,NumJ);                                                         % Wdij;
Wo = ones(1,NumJ);                                                         % Woij;
WO = 1;                                                                    
no = 0.2;
nd = 0.2;
ni = 0.2;
a = 0.9;
Rk = zeros(1,NumJ);Rk0 = zeros(1,NumJ);                                    % R(k),R(k-1) = [0],  % P(k),P(k-1) = [0],  % Q(k),Q(k-1) = [0];                             
Qk = zeros(1,NumJ);Qk0 = zeros(1,NumJ);                                    % ( R(1),R(0);Q(1),Q(0);P(1),P(0) );
Pk = zeros(NumI,NumJ);Pk0 = zeros(NumI,NumJ);

S = Is * Wi;                                                               % Sij;            ( S1j )                                     
X = Sigmoid(S);                                                            % Xij = f(Sij);   ( X1j )              
DfS = Dsigmoid(S);                                                         % f'(Sij);        ( f'(S1j) )
Os0 = X(1,:) * Wo';                                                        % O(k);           ( O(1) )

Rk = X(1,:) + Rk0;                                                         % R(k), P(k), Q(k) ;       
Qk = DfS(1,:) .* Qk0;                                                      % ( R(1), P(1), Q(1) ) ; 
for t1 = 1:NumI
    for t2 = 1:NumJ
        Pk(t1,t2) = Is(1,t1) * DfS(t2);
    end
end

Dwo = no * (Os0 - Os(1)) * Rk;                                             % DWo,Dwd,Dwo                      
Dwd = nd * (Os0 - Os(1)) * Wo .* Qk;                                       % ( DWo(1),Dwd(1),Dwo(1)  )
Dwi = zeros(size(Wi));
for t1 = 1 : NumI
    Dwi(t1,:) = ni * (Os0 - Os(1)) * Wo .* Pk(t1,:);
end

Rk0 = Rk;                                                                  % 存儲R(k), P(k), Q(k)為R(k-1), Q(k-1), P(k-1);
Qk0 = Qk;                                                                  
Pk0 = Pk;
Wo0 = Wo;                                                                  % 存儲Wo(k),Wd(k), Wi(k)為Wo(k-1),Wd(k-1),Wi(k-1);
Wd0 = Wd;
Wi0 = Wi;

Wo = Wo + Dwo;                                                             % 求W(k+1) ( W(2) );                                                            
Wd = Wd + Dwd;
Wi = Wi + Dwi;

for k = 2:size(Is,1)                                                       % 從第二個樣本的特征開始,到最后一個樣本的特征,循環
    
    S = Is * Wi;                                                           % Sij
    X = Sigmoid(S);                                                        % Xij = f(Sij) 
    DfS = Dsigmoid(S);                                                     % f'(Sij);
    Os0 = X(k,:) * Wo';                                                    % O(k)    

    error = Os0 - Os(k);                                                   % E(k)
    
    Rk = X(k,:) + Rk0;                                                     % R(k), P(k), Q(k) ;
    Qk = DfS(k,:) .* (X(k-1) + Qk0);
    for t1 = 1:NumI
        for t2 = 1:NumJ
            Pk(t1,t2) = Is(k,t1) * DfS(t2);
        end
    end
    Pk = Pk + Pk0;
     
    Dwo = no * (Os0 - Os(k)) * Rk;                                         % DWo,Dwd,Dwo   
    Dwd = nd * (Os0 - Os(k)) * Wo .* Qk;
    Dwi = zeros(size(Wi));
    for t1 = 1 : NumI
        Dwi(t1,:) = ni * (Os0 - Os(k)) * Wo .* Pk(t1,:);
    end

    Rk0 = Rk;                                                              % 存儲R(k), P(k), Q(k)為R(k-1), Q(k-1), P(k-1);
    Qk0 = Qk;
    Pk0 = Pk;
    Wo0 = Wo;                                                              % 存儲Wo(k),Wd(k), Wi(k)為Wo(k-1),Wd(k-1),Wi(k-1) ;
    Wd0 = Wd;
    Wi0 = Wi;
    
    Wo = Wo + Dwo;                                                         % 求W(k+1);
    Wd = Wd + Dwd;
    Wi = Wi + Dwi;
 
end

output = {Wo,Wd,Wi};