% alph, bet, b
% alph: M by [L by T] (cell)
% bet: M by [L by T] (cell)
% b: M by [L by T] (cell)
% P_M: M by 1
% Aij: L by L
% n: M by [N by T]
M = size(alph,1);
L = size(alph{1}(:,1),1);
pi_bar = 0;

P_M =
alph =
bet =
b =
Aij =


for m = 1:M
P_m = P_M(m);

% Calculate pi bar i
pi_bar = pi_bar + ( (1/P_m) * ( alph{m,1}(:,1) * bet{m,1}(:,1) ) );
end

% Calculate A bar ij
Aij = zeros(L,L);

% Calculate Lambda

% 10 neurons
lambda = zeros(L,10);

for j = 1:L
    for i = 1:L
        A_top = 0;
        A_bot = 0;
        for m = 1:M
            P_m = P_M(m);
            T_m = size(alph{m,1},2); % size of T
            for t = 1:(T_m - 1)
                A_top = A_top + (1/P_m) * alph{m,1}(i,t) * Aij(i,j) * b{m,1}(j,t+1) * bet{m,1}(j,t+1);
                A_bot = A_bot + (1/P_m) * alph{m,1}(j,t) * bet{m,1}(j,t+1);
            end
        end
        Aij(i,j) = A_top / A_bottom;
    end


    % lambda
    for k = 1:10
        lam_top = 0;
        lam_bot = 0;
        for m = 1:M
            P_m = P_M(m);
            T_m = size(alph{m,1},2); % size of T
            for t = 1:(T_m - 1)
                lam_top = lam_top + (1/P_m) * alph{m,1}(i,t) * n{m}(k,t+1) * bet{m,1}(j,t+1);
                lam_bot = lam_bot + (1/P_m) * alph{m,1}(j,t) * bet{m,1}(j,t+1);
            end
        end
        lambda(j,k) = lam_top / lam_bottom;
    end

end