Contents
function [Proba, Stat, Intern] = estimator_A(threshold,dim,varargin)
Parsing the arguments
error(nargoutchk(0,3,nargout));
p = inputParser;
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'}));
p.StructExpand = true;
p.KeepUnmatched = true;
p.parse(threshold,dim,varargin{:});
h = p.Results;
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);
N_rep = h.n-h.k+h.k_extra;
if nargout >1
Stat(length(h.threshold),1).rel_var = [];
Stat(length(h.threshold),1).rel_bias = [];
Stat(length(h.threshold),1).interval = [];
Stat(length(h.threshold),1).NbIter = [];;
Stat = orderfields(Stat);
end
if nargout==3
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);
dx = h.SCORE(X);
N = 1;
trial = h.n;
accepted = 0;
[thres_tmp Ind_vec] = HIGHER_SCORE(dx,h.k);
if nargout==3
thres_vec(1) = thres_tmp;
trial_vec(1) = h.n;
accept_vec(1) = 0;
prob_vec(1) = h.k/h.n;
mu_vec(1) = h.mu;
end
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;
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
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
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;
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);
Proba = k_prime/h.n*(h.k/h.n)^(N-1);
if nargout > 1
Stat = stats_estimator_A(h.n,h.k,'iter',N,'k_prime',k_prime);
end
else
Proba = 0;
if nargout > 1
Stat.rel_bias = 0;
Stat.rel_var = 0;
Stat.interval = [0,0];
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);