function y = mymvnpdfmultimulti(x, mu, sigma, doLog)
    if nargin < 4,
        doLog = false;
    end
    d = size(x);
    
    if length(d) == 4,
        x0 = bsxfun(@minus,x,mu);
        R = zeros([d(2),d(2),d(3),d(4)]);
        xRInv = zeros(d);
        for iPage = 1:d(3),
            for jPage = 1:d(4),
                R(:, :, iPage, jPage) = chol(sigma(:, :, iPage, jPage));
                xRInv(:, :, iPage, jPage) = x0(:, :, iPage, jPage) / R(:, :, iPage, jPage);
            end
        end
        logSqrtDetSigma = sum(log(R([1:4:((d(3)*d(4)-1)*4+1); 4:4:(d(3)*d(4)*4)])), 1);
        quadform = sum(xRInv.^2, 2);
%         y = exp(bsxfun(@minus, -0.5*permute(quadform, [3, 1, 4, 2]), reshape(logSqrtDetSigma, [d(3), 1, d(4)])) - d(2)*log(2*pi)/2);
        if doLog,
            y = bsxfun(@minus, -0.5*permute(quadform, [1, 3, 4, 2]), reshape(logSqrtDetSigma, [1, d(3), d(4)])) - d(2)*log(2*pi)/2;
        else
            y = exp(bsxfun(@minus, -0.5*permute(quadform, [1, 3, 4, 2]), reshape(logSqrtDetSigma, [1, d(3), d(4)])) - d(2)*log(2*pi)/2);
        end
    elseif length(d) == 3,
        x0 = bsxfun(@minus,x,mu);
        R = zeros([d(2),d(2),d(3)]);
        xRInv = zeros(d);
        for iPage = 1:d(3),
            R(:, :, iPage) = chol(sigma(:, :, iPage));
            xRInv(:, :, iPage) = x0(:, :, iPage) / R(:, :, iPage);
        end
        logSqrtDetSigma = sum(log(R([1:4:((d(3)-1)*4+1); 4:4:(d(3)*4)])), 1);
        quadForm = sum(xRInv.^2, 2);
%         y = exp(bsxfun(@minus, -0.5*permute(quadForm, [3, 1, 2]), logSqrtDetSigma') - d(2)*log(2*pi)/2);
        if doLog,
            y = bsxfun(@minus, -0.5*permute(quadForm, [1, 3, 2]), logSqrtDetSigma) - d(2)*log(2*pi)/2;
        else
            y = exp(bsxfun(@minus, -0.5*permute(quadForm, [1, 3, 2]), logSqrtDetSigma) - d(2)*log(2*pi)/2);
        end
    elseif length(d) == 2,
        d = size(x,2);
        x0 = bsxfun(@minus,x,mu);
        R = chol(sigma);
        xRInv = x0 / R;
        logSqrtDetSigma = sum(log(diag(R)));
        quadForm = sum(xRInv.^2, 2);
%         y = exp(-0.5*quadForm' - logSqrtDetSigma - d*log(2*pi)/2);
        if doLog,
            y = -0.5*quadForm - logSqrtDetSigma - d*log(2*pi)/2;
        else
            y = exp(-0.5*quadForm - logSqrtDetSigma - d*log(2*pi)/2);
        end
    end
end