getProbAggr.m

function [Pnir,logPnir] = getProbAggr(self, ist)
% [PNIR,LOGPNIR] = GETPROBAGGR(IST)
%
% Probability of a SPECIFIC per-class job distribution at a station.
% Returns P(n1 jobs of class 1, n2 jobs of class 2, ...) for current state.
%
% Compare with getProbMarg: returns queue-length distribution for a
% single class, i.e., P(n jobs of class r) for n=0,1,...,N(r).
%
% Input:
%   ist - Station index
%
% Output:
%   Pnir    - Scalar probability in [0,1]
%   logPnir - Log probability for numerical stability
 
if nargin<2 %~exist('ist','var')
    line_error(mfilename,'getProbAggr requires to pass a parameter the station of interest.');
end
sn = self.getStruct;
if ist > sn.nstations
    line_error(mfilename,'Station number exceeds the number of stations in the model.');
end
if isempty(self.result)
    self.run;
end
Q = self.result.Avg.Q;
N = sn.njobs;
if all(isfinite(N))
    switch self.options.method
        case 'exact'
            line_error(mfilename,'Exact marginal state probabilities not available yet in SolverMVA.');
        otherwise
            state = sn.state{sn.stationToStateful(ist)};
            [~, nir, ~, ~] = State.toMarginal(sn, ist, state);
            % Binomial approximation with mean fitted to queue-lengths.
            % Rainer Schmidt, "An approximate MVA ...", PEVA 29:245-254, 1997.
            logPnir = 0;
            for r=1:size(nir,2)
                logPnir = logPnir + nchoosekln(N(r),nir(r));
                logPnir = logPnir + nir(r)*log(Q(ist,r)/N(r));
                logPnir = logPnir + (N(r)-nir(r))*log(1-Q(ist,r)/N(r));
            end
            Pnir = real(exp(logPnir));
    end
else
    % Mixed or open model: use product-form distribution
    U = self.result.Avg.U;
    state = sn.state{sn.stationToStateful(ist)};
    [~, nir, ~, ~] = State.toMarginal(sn, ist, state);
    openClasses = find(isinf(N));
    closedClasses = find(isfinite(N));
    logPnir = 0;
 
    % Product-form probability for open classes
    if ~isempty(openClasses)
        if sn.sched(ist) == SchedStrategy.INF
            % Delay (infinite server): independent Poisson per class
            for r = openClasses
                if Q(ist,r) > 0
                    logPnir = logPnir + nir(r)*log(Q(ist,r)) - Q(ist,r) - gammaln(nir(r)+1);
                elseif nir(r) > 0
                    logPnir = -Inf;
                end
            end
        elseif sn.sched(ist) ~= SchedStrategy.EXT
            % Queue station: multinomial-geometric product form
            % P(n_1,...,n_R) = (1-rho) * n!/prod(n_r!) * prod(rho_r^n_r)
            rho_total = sum(U(ist, openClasses));
            n_total = sum(nir(openClasses));
            if rho_total < 1
                logPnir = logPnir + log(1 - rho_total) + gammaln(n_total + 1);
                for r = openClasses
                    rho_r = U(ist, r);
                    if nir(r) > 0
                        if rho_r > 0
                            logPnir = logPnir + nir(r)*log(rho_r) - gammaln(nir(r)+1);
                        else
                            logPnir = -Inf;
                        end
                    end
                end
            else
                logPnir = -Inf;
            end
        end
    end
 
    % Binomial approximation for closed classes
    % Rainer Schmidt, "An approximate MVA ...", PEVA 29:245-254, 1997.
    for r = closedClasses
        logPnir = logPnir + nchoosekln(N(r),nir(r));
        logPnir = logPnir + nir(r)*log(Q(ist,r)/N(r));
        logPnir = logPnir + (N(r)-nir(r))*log(1-Q(ist,r)/N(r));
    end
 
    Pnir = real(exp(logPnir));
end
end