Contents

function stat = stats_estimator_A(n,k,varargin)
% stat = stats_estimator_A(n,k,p)
% Gives some statistical properties of the estimator_A
%
% Inputs
% * n : total number of particles
% * k : number of survivors
%
% Parameters
% * alpha : confidence interval at alpha
%
% * P : probability to be estimated
% OR
% * iter : number of iterations
% * k_prime : number of particles above the last threshold
%
% Outputs
% stat : a cell containing
% * rel_bias : relative bias
% * rel_var : relative variance
% * interval : Interval of confidence
% * NbIter : number of iterations
%
% Comments
% -
%
% Authors
% * Implementation: Teddy Furon
% * Science: Frederic Cerou, Pierre Del Moral, Arnaud Guyader, Teddy Furon
%
% 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,k,varargin{:});
defaulted = p.UsingDefaults;
h = p.Results;
Input argument "n" is undefined.

Error in ==> stats_estimator_A at 57
p.parse(n,k,varargin{:});

Which mode?

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

Operational mode

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

    p0 = h.k/h.n;
    n0 = h.iter-1;
    r0 = h.k_prime/h.n;
    Prob = p0^n0*r0;
else

Theoretical mode

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

    p0 = h.k/h.n;
    n0 = floor(log(h.P)/log(p0));
    r0 = h.P*p0^(-n0);
    Prob = h.P;
end

Statistical properties

stat.NbIter = n0+1;
% relative bias of the estimation
stat.rel_bias = n0*(1-p0)/p0/h.n;
% relative variance of the estimation
stat.rel_var = (n0*(1-p0)/p0 + (1-r0)/r0)/h.n;
% Interval of alpha% confidence
La = sqrt(2)*erfinv(h.alpha);
stat.interval = Prob*(1-stat.rel_bias+sqrt(stat.rel_var)*[-La,La]);
stat.NbIter = n0+1;