Contents

function [Proba, Stat, Intern] = estimator_B(threshold,dim,varargin)
%[Proba, Stat, Intern] = estimator_B(threshold,dim,options)
%
% 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]
% * rate: acceptance rate                  [0.2]
% * verbose                                [true]
% * t: number of micro-replications        [20]
% * mu: strength of the micro-replication  [2]
% * mu_decay: decay rate of mu             [0.05]
% * max_iter: maximum number of iteration  [40000]
% * a: geometric weight of the average     [0.95]
% * N_sample: Sample rate of intern values
%
% Outputs
%
% * Proba: Estimated Probability
% * Stat : a structure containing some statistics about the estimation
% * Intern: Internal values for some iterations
%
% Comments
%
% * The replication process is done by t micro replications
% * This estimator is equivalent of estimator_A with k= n-1.
%
% Authors
%
% * Implementation: Teddy Furon
% * Science: Arnaud Guyader, Nicolas Hengartner, Eric Matzner-Lober
%
% Contact: teddy.furon@inria.fr
%
% Copyright INRIA 2010-2011
%
% 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)validateattributes(x,{'numeric'},...
    {'real', 'nonempty','finite'}));
p.addRequired('dim', @(x)validateattributes(x,{'numeric'},...
    {'integer', 'positive','scalar'}));

% Options
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.1,@(x)validateattributes(x,{'numeric'},...
    {'positive','scalar'}));
p.addParamValue('t',10,@(x)validateattributes(x,{'numeric'},...
    {'positive','scalar','integer'}));
p.addParamValue('max_iter',50000,@(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'},...
    {'positive','scalar'}));
p.addParamValue('verbose',true,@(x)islogical(x));
p.addParamValue('n',100*dim,@(x)validateattributes(x,{'numeric'},...
    {'scalar','integer', 'positive'}));
p.addParamValue('a',0.95,@(x)validateattributes(x,{'numeric'},...
    {'scalar', 'positive'}));
p.addParamValue('N_sample',1000,@(x)validateattributes(x,{'numeric'},...
    {'scalar','integer','positive'}));
p.addParamValue('minheap',false,@(x)islogical(x));

% 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 of estimator_B'
    disp(p.Results)
end
Input argument "dim" is undefined.

Error in ==> estimator_B at 76
p.addParamValue('n',100*dim,@(x)validateattributes(x,{'numeric'},...

Declare some variables

global dx indx; % This is required for inplace min_heap management.
threshold = max(h.threshold);
[target_thres_vec, I] = sort(h.threshold);
Proba  = zeros(size(h.threshold));
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 95% confidence
    Stat(length(h.threshold),1).NbIter = [];; % Nb of iterations
    Stat = orderfields(Stat);
end
if nargout==3
    % if Internal values are required
    Max_ind = ceil(h.max_iter /h.N_sample);
    thres_vec = zeros(Max_ind,1);
    trial_vec = zeros(Max_ind,1);
    accept_vec = zeros(Max_ind,1);
    prob_vec = zeros(Max_ind,1);
    mu_vec = zeros(Max_ind,1);
    iter_vec = zeros(Max_ind,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
rho = 1-1/h.n;
trial = h.n; % count number of detection call
thres_cnt = 1; % index of the next targeted threshold


% Init of the min-heap
if h.minheap
    indx = 1:h.n; % label these vectors
    min_heap; % sort dx (and indx) to become a min_heap.
    % the min is always the first element
    thres_tmp = dx(1); % define the first intermediate threshold
else
    [thres_tmp,Ind_min]=min(dx);
    Ind_rep = Ind_min;
end

% Init of the sliding average
sum_accepted = 0;
sum_a = 0;
iter_a = 0;

if nargout==3 % if Internal values are required
    thres_vec(1)=thres_tmp;
    trial_vec(1) = h.n; % count the number of detection trials
    accept_vec(1) = 0;
    prob_vec(1) = rho;
    mu_vec(1) = h.mu;
    iter_vec(1) = 1;
    N_internal = 1;
end

% Display
if h.verbose
    fprintf('Iteration \tThreshold \tProb \t\tmu \t\tTrials \n');
    fprintf('%d\t\t%6.2e\t%6.2e\t%6.2e\t%6.2e\n',[N,thres_tmp,rho.^(N-1),h.mu,trial]);
end

Iterate

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

Draw a particule to be replicated

    if h.minheap
        Ind_heap = 1+ceil((h.n-1)*rand(1,1));% anyone but the first (min)
        Ind_rep = indx(Ind_heap);
        dz_a = dx(Ind_heap);
    else
        while Ind_rep==Ind_min
            Ind_rep = ceil(h.n*rand(1));
        end
        dz_a = dx(Ind_rep);
    end
    z_a = X(:,Ind_rep);

    accepted = 0;

Replication process

    for iT = 1:h.t % h.t consecutive micro-replications
        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
    if h.minheap
        X(:,indx(1))= z_a;
        min_heap(dz_a); % integrate the new score in the min_heap
        thres_tmp  = dx(1); % new minimum
    else
        X(:,Ind_min) = z_a;
        dx(Ind_min) = dz_a;
        [thres_tmp,Ind_min]=min(dx);
        Ind_rep = Ind_min;
    end

is the strength too big?

    sum_accepted = accepted + h.a*sum_accepted; %  sum over sliding wind.
    sum_a = sum_a+h.a^iter_a; % sum of the weights
    iter_a = iter_a +1; % length of the sliding window

    if ((sum_accepted / sum_a)<=(h.rate*h.t))&(iter_a>10)
        h.mu = h.mu*(1-h.mu_decay); % the strength of the replication is too big
        sum_accepted = 0; % the sliding average is reset.
        sum_a = 0;
        iter_a = 0;
    end

Have we reached the next threshold?

    if thres_tmp>target_thres_vec(thres_cnt)
        % Estimated probability
        Proba(I(thres_cnt)) = rho^(N-1);
        if nargout>1
            % Statistics
            Stat(I(thres_cnt)) = stats_estimator_B(h.n,'iter',N);
        end
        thres_cnt = thres_cnt + 1;
    end

Internal values

    if (mod(N,h.N_sample)==0)

Keep records for statistics

        if nargout ==3
            N_internal  = N_internal + 1;
            iter_vec(N_internal) = N;
            thres_vec(N_internal) = thres_tmp;
            prob_vec(N_internal) = rho^N;
            mu_vec(N_internal) = h.mu;
            trial_vec(N_internal) = trial;
        end

Display some data

        if (h.verbose)
            fprintf('%d\t\t%6.2e\t%6.2e\t%6.2e\t%6.2e\n',[N,thres_tmp,rho.^(N-1),h.mu,trial]);
        end
    end
end

Conclude

if N<h.max_iter
    % Estimated probability
    Proba(I(end)) = rho^(N-1);
    % Statistics
    if nargout>1
        Stat(I(end)) = stats_estimator_B(h.n,'iter',N);
    end
else % failure, too few iterations
    % Estimated probability
    Proba(I(end)) = 0;
    % Statistics
    if nargout >1
        Stat(I(end)).NbIter = h.max_iter;
    end
end

Output internal variables

if nargout==3
    Intern = [iter_vec(1:N_internal), thres_vec(1:N_internal),prob_vec(1:N_internal),trial_vec(1:N_internal),mu_vec(1:N_internal)];
    Intern = [Intern;N,thres_tmp,Proba(I(end)),trial,h.mu];
end

SubFunctions - min heap

min_heap, min_heap_el

function min_heap(x)

Min_heap

Maintain the heap so that the min is always the first element

Global variables * dx is a vector of scores * indx is a vector of indices

No input * Build the min-heap

One input * x : a new score to be inserted in the min-heap

global dx indx

Init - build the min-heap

if nargin==0
    % there is a faster way but a sort will do the job.
    [dx,indx] = sort(dx,'ascend');
else

Insert a new data

assume that dx is already a min-heap the job is now to insert x

    dx(1) = x;
    min_heap_el(1); % recursive elementary operation
end
function min_heap_el(i)

Elementary operation to maintain the min_heap

from wikipedia http://en.wikipedia.org/wiki/Binary_heap

global dx indx
ind_left = 2*i;
ind_right = 2*i+1;

if (ind_left <= length(dx))&(dx(ind_left)<dx(i))
    ind_smallest = ind_left;
else
    ind_smallest = i;
end
if (ind_right<=length(dx))&(dx(ind_right)<dx(ind_smallest))
    ind_smallest=ind_right;
end
if ind_smallest~=i
    % swap
    dx_tmp = dx(ind_smallest);
    ind_tmp = indx(ind_smallest);

    dx(ind_smallest) = dx(i);
    indx(ind_smallest) = indx(i);

    dx(i) = dx_tmp;
    indx(i) = ind_tmp;

    min_heap_el(ind_smallest);
end

SubFunctions - Default problem

SCORE_INTERNAL, GENERATE_INTERNAL, MODIFY_INTERNAL
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);