...

function [mean,std,muLN,sigmaLN] = ...
    PoissRatioStat(lambdaNum,lambdaDen,constant)
% Mean and standard deviation 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. 
%
% Syntax:
% [m,v] = PoissRatioStat(x,lambdaNum,lambdaDen)
% or
% [m,std,muLNfit,sigmaLNfit] = PoissRatioStat(lambdaNum,lambdaDen,constant)
%
% Example (Compute the exact mean and variance of X and the parameters
% of the fitted Log-normal distribution):
%  lambdaNum = 10;
%  lambdaDen = 2000;
%  constant = 300;
%  [m,s,mLN,sLN] = PoissRatioStat(lambdaNum,lambdaDen,constant)
%
% See also: PoissRatioCdf, PoissRatioPmf
%
% 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 < 3
    constant = [];
end

if nargin < 2 || nargin > 3,
    error('VW:PoissRatioStat:TooFewInputgas',...
        'Requires two input arguments.');
end

if ~isreal(lambdaNum) || ~isreal(lambdaDen) || lambdaNum <= 0 || ...
        lambdaDen <= 0 || ~isscalar(lambdaNum) || ~isscalar(lambdaDen)
    error('VW:PoissRatioStat:NegativeParameters',...
        'Parameters lambdaNum and lambdaDen must be real and positive scalars.');
end
if isempty(constant)
    constant = 1;
end
%% Expectation and standard deviation
e = (1- exp(-lambdaDen)) * lambdaNum / lambdaDen;
s = stdN(lambdaNum,lambdaDen);
mean = constant * e;
std = constant * s;
%% Parameters of the fitted lognormal distribution
sigmaLN = sqrt(log(1 + (std / mean)^2));
muLN = log(mean) - sigmaLN^2/2;
%% function stdN / regularized standard deviation
%Computes the exact regularized standard deviation of the ratio of
%independent Poisson random variables

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

function s = stdN(lambdaNum,lambdaDen)
eulerGamma = 5.772156649015329e-001;
ratio = lambdaNum / lambdaDen;
aux1 = - ((1 - exp(-lambdaDen)) * ratio)^2;
aux2 = - (1 + lambdaNum ) * ratio;
aux3 = aux2 * exp(-lambdaDen) * (eulerGamma + log(lambdaDen));
if lambdaDen < 700,
    aux4 = aux2 * real(-expint(-lambdaDen)) * exp(-lambdaDen);
else
    aux4 = aux2 / (lambdaDen - 1);
end
s = sqrt(aux1 + aux3 - aux4);