Contents

function [q, Stat, Intern] = estimator_quantile_B_cell(P,dim,varargin)
%[q, Stat, Intern] = estimator_quantile_B_cell(P,dim,varargin)
%
% /!\ This implementation uses cells!
%
% Inputs
% * P : probability
% * 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             [1.1]
% * max_iter: maximum number of iteration  [40000]
% * a: geometric weight of the average     [0.95]
% * N_sample: Sample rate of intern values [1000]
% * alpha: a 100*alpha% confidence interval[0.95]
%
% 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
% * This estimator is equivalent of estimator_A with k= n-1.
% * Works with cells
%
% 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('P', @(x)validateattributes(x,{'numeric'},...
    {'real', 'nonempty','finite','positive'}));
p.addRequired('dim', @(x)validateattributes(x,{'numeric'},...
    {'integer', 'positive'}));

% 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',40000,@(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('alpha',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(P,dim,varargin{:});

h = p.Results;

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

Error in ==> estimator_quantile_B_cell 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

La = sqrt(2)*erfinv(h.alpha); % quantile for alpha% confidence interval
q = zeros(size(P)); % estimated quantiles
P = P(:);
m_tmp = repmat(1:3,length(P),1);
m_mat = [zeros(3*length(P),1),repmat((1:length(P))',3,1),m_tmp(:)];
m_mat(1:length(P),1) = floor(-h.n*log(P)-La*sqrt(-h.n*log(P))); %Lower bound
m_mat((1:length(P))+length(P),1) = ceil(log(P)/log(1-1/h.n));  %Estimate
m_mat((1:length(P))+2*length(P),1)=ceil(-h.n*log(P)+La*sqrt(-h.n*log(P)));%Upper bound
% m_mat stores:
% 1st column: the index where to stop to estimate the corresponding
% quantiles
% 2nd column: the index of the corresponding probability in the input
% vector
% 3rd column: the nature of the quantile (Lower bound, Estimates, Upper)

[target_m_vec,I] = sort(m_mat(:,1));
m_p = target_m_vec(end); % Biggest level
if nargout > 1
    Stat(length(P),1).NbIter = [];
    Stat(length(P),1).interval = [];
    Stat = orderfields(Stat);
end
if nargout==3 % if Internal values are required
    Max_ind = ceil(m_p /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); % cell vector 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 calls
m_cnt = 1; % counter for targeted level

% 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

for N = 2:m_p
    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

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

Have we reach the next targeted levels ?

    if N==target_m_vec(m_cnt)
        m_ind = find(N==target_m_vec); % there might be more than 1 match
        ind_cnt = I(m_ind);
        for ki=1:length(ind_cnt)
            k = ind_cnt(ki);
            if m_mat(k,3)==2 % EStimate
                q(m_mat(k,2)) = thres_tmp;
            end
            if nargout >1
                switch m_mat(k,3)
                    case 1 % Lower bound of an IC
                        Stat(m_mat(k,2)).interval(1)=thres_tmp;
                        Stat(m_mat(k,2)).NbIter(1)=N;
                    case 2 % Estimate
                        Stat(m_mat(k,2)).NbIter(2)=N;
                    case 3 % Upper bnd of an IC
                        Stat(m_mat(k,2)).interval(2)=thres_tmp;
                        Stat(m_mat(k,2)).NbIter(3)=N;
                end
            end
        end
        m_cnt = m_cnt+length(m_ind); % Go to the next levels
    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,(1-1/h.n)^N,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
[n e] = size(X);
d = zeros(n,1);
for k=1:n
    d(k) = abs(X{k}(1)./sqrt(sum([X{k}].^2,1)));
end
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 = cell(n,1);
for k=1:n
    X{k} = randn(dim,1);
end
function y = MODIFY_INTERNAL(x,mu)

MODIFY_INTERNAL

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