Contents
- Parsing the arguments
- Declare some variables
- Initial step
- Iterate
- Draw a particule to be replicated
- Replication process
- is the strength too big?
- Internal values
- Keep records for statistics
- Display some data
- Have we reach the next targeted levels ?
- Output internal variables
- SubFunctions - min heap
- Min_heap
- Init - build the min-heap
- Insert a new data
- Elementary operation to maintain the min_heap
- SubFunctions - Default problem
- SCORE_INTERNAL
- GENERATE_INTERNAL
- MODIFY_INTERNAL
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);