Умный и быстрый способ вычислить тройное суммирование?

Всем привет,

Мне нужно вычислить следующую тройную сумму для разных значений n для всех s.

Также известны ~ p, l_s и q_s1, q_s2 и q_s3, а π — полиномиальное распределение. Это довольно легко сделать с помощью циклов for (что я сделал ниже), но выполнение занимает вечность по мере увеличения n. Для сравнения, при n=10 это занимает 345 секунд. Мне нужно работать с n=10000.

Есть ли умный способ ускорить процесс?

clear
clc
load data.mat
tic
n = 10;
l_s = height( segments );       % # of segments
L = sum( segments.NewLength );      % total legnth of Segments
expSlength = zeros(l_s,1);
meanPass = median( barriers.prePass );
qs1 = segments.NewLength ./ L;
qs2 = segments.segDis2Mth ./ L;
qs3 = (L - segments.NewLength - segments.segDis2Mth) ./ L;
for s = 1 : l_s
    for k = 0 : n
        innersum = 0;
        for t = 0 : n-k
            B = 0;
            for r = 0 : k
                A = segments.NewLength(s) / (k + 1) * meanPass^r;
                B = B + A;
            end
            PROB = [qs1(s), qs2(s), qs3(s)];
            X = [ k, t, n-k-t ];
            p_stk = mnpdf(X, PROB);
            innersum = innersum + p_stk * meanPass^t * B;
        end
        expSlength(s) = expSlength(s) + innersum;
    end
end
toc

ПРИМЕЧАНИЕ. 

Matlabsolutions.com предоставляет последнюю Помощь по домашним заданиям MatLab, Помощь по заданию MatLab для студентов, инженеров и исследователей в различных отраслях, таких как ECE, EEE, CSE, Mechanical, Civil со 100% выходом. Код Matlab для BE, B.Tech ,ME,M.Tech, к.т.н. Ученые со 100% конфиденциальностью гарантированы. Получите проекты MATLAB с исходным кодом для обучения и исследований.

ЧАСТЬ 4: Петли заменены местами:

tic
n   = 50; 
l_s = numel(segments.batNetID);   % # of segments
L   = sum(segments.NewLength);    % total legnth of Segments
expSlength = zeros(l_s, 1);
meanPass = median( barriers.prePass );
qs1 = segments.NewLength ./ L;
qs2 = segments.segDis2Mth ./ L;
qs3 = (L - segments.NewLength - segments.segDis2Mth) ./ L;
meanPass_pow  = cumprod([1, repmat(meanPass, 1, n)]);
meanPass_powC = cumsum(meanPass_pow);
gammaln2 = gammaln(1:n+1);
logPROB  = log([qs1, qs2, qs3]).';
for k = 0:n
   X = zeros(n - k + 1, 3);
   for t = 0:n-k
      X(t+1, :) = [k, t, n-k-t];
   end
   c = gammaln2(n + 1) - sum(gammaln2(X + 1), 2);
     % VERSION 1: full inner loop
     % for s = 1:l_s
     %   B = meanPass_powC(k + 1) * segments.NewLength(s) / (k + 1);      
     %   p_stk = exp(c + X * logPROB(:, s));
     %   
     %   innersum = 0;
     %   for t = 0 : n-k
     %      innersum = innersum + p_stk(t+1) * meanPass_pow(t+1);
     %   end
     %   expSlength(s) = expSlength(s) + innersum * B;
     % end
     % VERSION 2: Partially vectorized
     % p_stk = exp(c + X * logPROB);
     % B = segments.NewLength * (meanPass_powC(k + 1) / (k + 1));
     % for s = 1:l_s
     %    innersum = 0;
     %    for t = 0 : n-k
     %       innersum = innersum + p_stk(t+1, s) * meanPass_pow(t+1);
     %    end
     %    expSlength(s) = expSlength(s) + innersum * B(s);
     % end
     % VERSION 3: fully vectorized
     p_stk      = exp(c + X * logPROB);
     B          = segments.NewLength * (meanPass_powC(k + 1) / (k + 1));
     expSlength = expSlength + (meanPass_pow(1:n-k+1) * p_stk).' .* B;  
  end
  toc
 n = 50;
 ARRAYFUN:   8.2 sec
 VERSION 1: 19.5 sec
 VERSION 2:  6.2 sec
 VERSION 3:  5.0 sec

СМОТРИТЕ ПОЛНЫЙ ОТВЕТ НАЖМИТЕ НА ССЫЛКУ