Contents

function [Proba, Stat, Intern] = estimator_A(threshold,dim,varargin)
% [Proba, Stat, Intern] = estimator_A(threshold,dim,option)
%
% Inputs
% * threshold : defines the acceptance region score(x)>threshold
% * dim : Dimension of the particule space
% * option: a structure containing the following parameters (fields)
%
% Parameters
% * GENERATE: particule generator          [@(x)GENERATE_INTERNAL]
% * SCORE: particule score                 [@(x)SCORE_INTERNAL]
% * MODIFY: particule modifier             [@(x)MODIFY_INTERNAL]
% * n: total number of particules          [100*dim]
% * k: number of survivors                 [50*dim]
% * rate: acceptance rate                  [0.7]
% * verbose                                [false]
% * t: number of micro-replications        [10]
% * mu: strength of the micro-replication  [10]
% * mu_decay: decay rate of mu             [0.05]
% * max_iter: maximum number of iteration  [100]
%
% Outputs
% * Proba: Estimated Probability
% * Stat : a structure containing some statistics about the estimation
% * Intern: Internal values per iteration
%
% Comments
% * The replication process is done by T micro replications
% * only the n-k+k_extra particles are replicated.
%
% 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('threshold', @(x)isscalar(x));
p.addRequired('dim', @(x)validateattributes(x,{'numeric'},{'integer', 'positive'}));

p.addParamValue('GENERATE',@GENERATE_INTERNAL,@(x)isa(x,'function_handle'));
p.addParamValue('SCORE',@SCORE_INTERNAL,@(x)isa(x,'function_handle'));
p.addParamValue('MODIFY',@MODIFY_INTERNAL,@(x)isa(x,'function_handle'));
p.addParamValue('rate',0.3,@(x)validateattributes(x,{'numeric'},{'nonnegative','scalar'}));
p.addParamValue('t',10,@(x)validateattributes(x,{'numeric'},{'positive','scalar','integer'}));
p.addParamValue('max_iter',200,@(x)validateattributes(x,{'numeric'},{'positive','scalar','integer'}));
p.addParamValue('mu',2,@(x)validateattributes(x,{'numeric'},{'positive','scalar'}));
p.addParamValue('mu_decay',0.05,@(x)validateattributes(x,{'numeric'},{'nonnegative','scalar'}));
p.addParamValue('verbose',true,@(x)islogical(x));
p.addParamValue('n',100*dim,@(x)validateattributes(x,{'numeric'},{'integer', 'positive'}));
p.addParamValue('k',50*dim,@(x)validateattributes(x,{'numeric'},{'integer', 'positive'}));
p.addParamValue('k_extra',0,@(x)validateattributes(x,{'numeric'},{'integer', 'positive'}));

% Parse and validate all input arguments.
p.StructExpand = true;
p.KeepUnmatched = true;
p.parse(threshold,dim,varargin{:});

h = p.Results;

% Display the names of all arguments.
if h.verbose
    fprintf('\n')
    disp 'List of all arguments:'
    disp(p.Results)
end
Input argument "dim" is undefined.

Error in ==> estimator_A at 62
p.addParamValue('n',100*dim,@(x)validateattributes(x,{'numeric'},{'integer', 'positive'}));

Declare some variables

Y = zeros(h.dim,h.k); % matrix of survivor vectors
N_rep = h.n-h.k+h.k_extra; % Number of replicated particles

if nargout >1
    % Create the structure
    Stat(length(h.threshold),1).rel_var = []; % relative variance of the estimation
    Stat(length(h.threshold),1).rel_bias = [];% relative bias of the estimation
    Stat(length(h.threshold),1).interval = [];% Interval of confidence
    Stat(length(h.threshold),1).NbIter = [];; % Nb of iterations
    Stat = orderfields(Stat);
end
if nargout==3 % if Internal values are required
    thres_vec = zeros(h.max_iter,1);
    trial_vec = zeros(h.max_iter,1);
    accept_vec= zeros(h.max_iter,1);
    prob_vec  = zeros(h.max_iter,1);
    mu_vec = zeros(h.max_iter,1);
end

Initial step

X = h.GENERATE(h.dim,h.n); % matrix of generated vectors
dx =  h.SCORE(X); % calculate their score
N = 1; % Number of iteration
trial = h.n; % Number of detection trials
accepted = 0;
[thres_tmp Ind_vec] = HIGHER_SCORE(dx,h.k); % take the k-th higher score

if nargout==3 % if Internal values are required
    thres_vec(1) = thres_tmp; % define the first intermediate threshold
    trial_vec(1) = h.n;
    accept_vec(1) = 0;
    prob_vec(1) = h.k/h.n;
    mu_vec(1) = h.mu;
end

% Display
if h.verbose
    fprintf('Iteration \tThreshold \tProb \t\tTrials \t\tAccept \t\tmu \n');
    fprintf('%6.2e\t%6.2e\t%6.2e\t%6.2e\t%6.2e\t%6.2e\n',[N thres_tmp,(h.k/h.n)^N,trial,accepted/N_rep/h.t,h.mu]);
end

Iterate

while(thres_tmp<h.threshold)&(N<h.max_iter)
    N = N+1;

select the survivors whose score > intermediate threshold

    Y=X(:,Ind_vec);
    dy = dx(Ind_vec);

Select the indices of the vectors to be replicated

    permut = randperm(h.k);
    Index_vec = permut(mod((1:N_rep),h.k)+1);

Replication process

    flag_valid_iter = false;
    while ~flag_valid_iter
        trial = trial + (N_rep)*h.t;
        accepted = 0; % counts the number of successful replications
        for i=1:N_rep
            Index = Index_vec(i);
            z_a = Y(:,Index);
            dz_a = dy(Index);
            flag_valid_micro = false;
            for iT = 1:h.t % micro-replication
                z_b = h.MODIFY(z_a,h.mu);
                dz_b = h.SCORE(z_b);
                if dz_b >=thres_tmp
                    z_a = z_b;
                    dz_a = dz_b;
                    accepted = accepted +1;
                end
            end % end micro-replication

            X(:,i)= z_a;
            dx(i) = dz_a;
        end
        Index_rem = permut(mod((N_rep+1:h.n),h.k)+1);
        X(:,(N_rep+1):h.n) = Y(:,Index_rem);
        dx((N_rep+1):h.n) = dy(Index_rem);
        if (accepted>(h.rate*N_rep*h.t))
            flag_valid_iter = true;
            % iteration is correct when rate > acceptance rate%
        else
            h.mu = h.mu*(1-h.mu_decay);
        end
    end
    [thres_tmp Ind_vec] = HIGHER_SCORE(dx,h.k);

Keep records for statistics

    if nargout==3
        accept_vec(N) = accepted/N_rep/h.t;
        mu_vec(N) = h.mu;
        thres_vec(N) = thres_tmp;
        prob_vec(N) = (h.k/h.n)^N;
        trial_vec(N) = trial;
    end

Display some data

    if (h.verbose)
        fprintf('%6.2e\t%6.2e\t%6.2e\t%6.2e\t%6.2e\t%6.2e\n',[N thres_tmp,(h.k/h.n)^N,trial,accepted/N_rep/h.t,h.mu]);
    end
end

Conclude

if N<h.max_iter
    Ind_final = find(dx>h.threshold);
    k_prime = length(Ind_final);
    % Estimated probability
    Proba = k_prime/h.n*(h.k/h.n)^(N-1);
    % Statistics
    if nargout > 1
        Stat = stats_estimator_A(h.n,h.k,'iter',N,'k_prime',k_prime);
    end
else % failure, too few iterations
    % Estimated probability
    Proba = 0;
    % Statistics
    if nargout > 1
        Stat.rel_bias = 0;      % relative bias of the estimation
        Stat.rel_var = 0;       % relative variance of the estimation
        Stat.interval = [0,0];  % Interval of confidence
        Stat.NbIter = h.max_iter;
    end
end

Output internal variables

if nargout==3
    N_vec = 1:N;
    Intern=[N_vec',thres_vec(N_vec),prob_vec(N_vec),...
        trial_vec(N_vec),accept_vec(N_vec),mu_vec(N_vec)];
end

Subfunctions

HIGHER_SCORE, SCORE_INTERNAL, GENERATE_INTERNAL, MODIFY_INTERNAL
function [t Ind_vec] = HIGHER_SCORE(dx,k)

HIGHER SCORE

[t Ind_vec] = HIGHER_SCORE(dx,k)
values of the k-th higher scores
Inputs
* dx : a vector of scores
* k : integer such that length(dx)>k
Outputs
* t : k-th higheest score
* Ind_vec : vector of k indices giving the kth highest scores in dx
if length(dx)<k
    error('k greater than the length of dx');
end

if k~=1
    [vec I] = sort(dx,'descend');
    t = vec(k);
    Ind_vec = I(1:k);
else


end
function d = SCORE_INTERNAL(X)

SCORE_INTERNAL

d = SCORE(X)
Calculate the score for all the vectors store in matrix X
the score is here the absolute value of the  normalized correlation
* X : a vector or a matrix
* d : a scalar or a vector
d = abs(X(1,:)./sqrt(sum(X.^2,1)));
function X = GENERATE_INTERNAL(dim,n)

GENERATE_INTERNAL

X = GENERATE(L,n)
Generate n vectors of length dim
distributed as p_X (here, white gaussian noise)
X = randn(dim,n);
function y = MODIFY_INTERNAL(x,mu)

MODIFY_INTERNAL

y = MODIFY_INTERNAL(x,mu)
replication process for white gaussian vector
* x : vector to be replicated
* mu : strength of the replication
* y : replicated vector
y = (x + mu*randn(size(x)))/sqrt(1+mu^2);