PoissRatioPmf(x,lambdaNum,lambdaDen,constant,...

function [pmf,x] = PoissRatioPmf(x,lambdaNum,lambdaDen,constant,...
    isPlot,critfactor)
% Probability mass function (pmf) of the random variable X,
% X = constant * (Y_lambdaNum / Y_lambdaDen), which is proportional to the
% ratio of two independent Poisson random variables: (i) Y_lambdaNum, with
% parameter lambdaNum, and (ii) Y_lambdaDen which follows the 1-displaced
% Poisson distribution (shifted by 1 to the right), with parameter
% lambdaDen, respectively. Then, the probability mass function  is
% pmf = Prob(constant * (Y_lambdaNum / Y_lambdaDen) = x).
%
% Syntax:
% pmf = PoissRatioPmf(x,lambdaNum,lambdaDen);
% or
% [pmf,x] = PoissRatioPmf([],lambdaNum,lambdaDen);
% or
% [pmf,x] = PoissRatioPmf(x,lambdaNum,lambdaDen,constant,isPlot,critfactor);
%
% Here, isPlot is logical flag for ploting the pmf (default value:
% isPlot = false).
% The parameter 'critfactor' sets the heuristic rule
% 'expectation(X) < critfactor * standardDeviation(X)' to choose
% either the exat method for computing the pmf, or the log-normal
% approximation (default value: critfactor = 10). If 'critfactor = Inf' the
% algorithm calculates the cdf by the exact method (this could take long
% time for large values of lambdaDen / even if the simple log-normal
% approximation is almost perfect).
%
% Example (Compute the exact pmf in automatically selected values of x):
%  lambdaNum = 2;
%  lambdaDen = 5;
%  constant = 1;
%  isPlot = true;
%  critfactor = 10;
%  [pmf,x] = PoissRatioPmf([],lambdaNum,lambdaDen,constant, ...
%     isPlot,critfactor);
%
% See also: PoissRatioCdf
%
% References:
% ARENDACK, B. - SCHWARZ, K. - TOLC JR, S. - WIMMER, G. - WITKOVSK, V.:
% Variability issues in determining concentration of isoprene in human
% breath by PTR-MS. Journal of Breath Research 2, 2008, 037007 (8pp),
% doi:10.1088/1752-7155/2/3/037007

% (c) Viktor Witkovsky ([e-mail])
% Revised: 14-Nov-2009 15:36:29
%% Check the inputs
if nargin < 6
    critfactor = [];
end

if nargin < 5
    isPlot = [];
end

if nargin < 4
    constant = [];
end

if nargin < 3 || nargin > 6,
    error('VW:PoissRatioPmf:TooFewInputgas',...
        'Requires three input arguments.');
end

if ~isreal(lambdaNum) || ~isreal(lambdaDen) || lambdaNum <= 0 || ...
        lambdaDen <= 0 || ~isscalar(lambdaNum) || ~isscalar(lambdaDen)
    error('VW:PoissRatioPmf:NegativeParameters',...
        'Parameters lambdaNum and lambdaDen must be real and positive scalars.');
end
if isempty(constant)
    constant = 1;
end

if isempty(isPlot)
    isPlot = false;
end

if isempty(critfactor)
    critfactor = 10;
end
%% Calculate the pmf at x
if isempty(x)
    N = 1000;
    stdTol = 4;
    x = SetX(lambdaNum,lambdaDen,constant,stdTol,N);
end

epsilon = 100 * eps;
pmf = PoissRatioCdf(x,lambdaNum,lambdaDen,constant,[],critfactor);
pmf = pmf - PoissRatioCdf(x - epsilon,lambdaNum,lambdaDen,...
    constant,[],critfactor);
%% Plot the pmf
if isPlot
    [X,ind] = unique(x);
    bar(X,pmf(ind))
    title('Exact PMF of the Ratio of Independent Poisson Random Variables')
    xlabel('x')
    ylabel('Probability Mass Function')
end
%% function SetX
function x = SetX(lambdaNum,lambdaDen,constant,stdTol,N)
% Set the best N rational numbers as possible outcomes of the ratio of
% two Poisson random variables with intensities lambdaNum and lambdaDen
% X = constant * Y_lambdaNum / Y_lambdaDen
%
% Syntax
% x = SetX(lambdaNum,lambdaDen,constant,N,stdTol)

% (c) Viktor Witkovsky ([e-mail])
% Revised: 14-Nov-2009 15:36:29

stdNum = ceil(sqrt(lambdaNum));
lowNum = max(0,lambdaNum - stdTol * stdNum);
uppNum = lambdaNum + stdTol * stdNum;
stdDen = ceil(sqrt(lambdaDen));
lowDen = max(1,lambdaDen - stdTol * stdDen);
uppDen = lambdaDen + stdTol * stdDen;
indNum = lowNum:uppNum;
indDen = lowDen:uppDen;
[x,i,j] = unique(kron(indNum',1./indDen));
if length(x) > N
    jSort=sort(j);
    [jSortUnique,iInd] = unique(jSort);
    [iInddiff,iSort] = sort(diff([0;iInd]),'descend');
    ind = iSort(1:N);
    x = sort(x(jSort(iInd(ind))));
end
[num,den] = rat(x,1.e-6);
x = unique(num ./ den);
x = constant * x;