Contents

function stat = stats_estimator_B(n,varargin)
% stat = stats_estimator_B(n,options)
% Gives some statistical properties of the estimator_B
%
% Inputs
% * n : number of particles
%
% Options
% * alpha : confidence interval at alpha
%
% * P : probability to be estimated
% OR
% * iter : number of iterations
%
% Outputs
% * stat : a struct containing
% * rel_bias : relative bias
% * rel_var : relative variance
% * interval : Interval of confidence
% * NbIter : number of iterations
%
% Comments
% -
% Authors
% * Implementation: Teddy Furon
% * Science: Arnaud Guyader, Nicolas Hengartner, Eric Matzner-Lober
%
% Contact: teddy.furon@inria.fr
%
% Copyright INRIA 2010
%
% Licence
% This software and its subprograms are governed by the CeCILL license
% under French law and abiding by the rules of distribution of free software.
% See http://www.cecill.info/licences.en.html

Parsing the arguments

error(nargoutchk(0,3,nargout));

% Create an instance of the inputParser class.
p = inputParser;
% Define inputs that one must pass on every call:
p.addRequired('n', @(x)validateattributes(x,{'numeric'},{'integer', 'positive'}));
%p.addRequired('k', @(x)validateattributes(x,{'numeric'},{'integer', 'positive'}));

%p.addParamValue('k_prime',1,@(x)validateattributes(x,{'numeric'},{'integer', 'positive'}));
p.addParamValue('iter',100,@(x)validateattributes(x,{'numeric'},{'integer', 'positive'}));
p.addParamValue('P',10^-7,@(x)validateattributes(x,{'numeric'},{'positive','scalar'}));
p.addParamValue('alpha',0.95,@(x)validateattributes(x,{'numeric'},{'positive','scalar'}));

% Parse and validate all input arguments.
p.StructExpand = true;
p.KeepUnmatched = true;
p.parse(n,varargin{:});
defaulted = p.UsingDefaults;
h = p.Results;
Input argument "n" is undefined.

Error in ==> stats_estimator_B at 55
p.parse(n,varargin{:});

Which mode?

if nnz(strcmp('P',defaulted))>0

Operational mode

'P' was not passed a real estimator was runned before

    n0 = h.iter-1;
    Prob = (1-1/h.n)^n0;
else

Theoretical mode

'P' was passed We want an estimate of the properties for such a probability to be estimated

    n0 = floor(log(h.P)/log(1-1/h.n));
    Prob = h.P;
end

Statistical properties

Number of iteration

stat.NbIter = n0+1;
% relative bias of the estimation
stat.rel_bias = 0;
% relative variance of the estimation
stat.rel_var = (Prob^(-1/h.n)-1);
% Interval of alpha% confidence
La = sqrt(2)*erfinv(h.alpha);
stat.interval = Prob*exp(La*[-1,1]'/sqrt(h.n)*sqrt(-log(Prob)+La^2/4/h.n)-...
    La^2/2/h.n);
% Order fields
stat = orderfields(stat);