Test for estimator_A

Estimation of a weak probability

Authors

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

Contents

Example

What is the probability that

Data

dim = 20;
threshold = 0.95;
n = 2000; % Total number of particles
k = 1000; % number of survivors

% real value for this problem
Proba_true = 1 - betainc(threshold^2,1/2,(dim-1)/2);

Prediction of the Estimator performance

Stat_pre = stats_estimator_A(n,k,'P',Proba_true)
Stat_pre = 

      NbIter: 35
    rel_bias: 0.0170
     rel_var: 0.0171
    interval: [3.4177e-11 5.8303e-11]

Estimation

[Proba_est, Stat, Intern] = estimator_A(threshold,dim,'n',n,'k',k);
List of all arguments:
          dim: 20
     GENERATE: @GENERATE_INTERNAL
            k: 1000
      k_extra: 0
     max_iter: 200
       MODIFY: @MODIFY_INTERNAL
           mu: 2
     mu_decay: 0.0500
            n: 2000
         rate: 0.3000
        SCORE: @SCORE_INTERNAL
            t: 10
    threshold: 0.9500
      verbose: 1

Iteration 	Threshold 	Prob 		Trials 		Accept 		mu 
1.00e+00	1.54e-01	5.00e-01	2.00e+03	0.00e+00	2.00e+00
2.00e+00	2.59e-01	2.50e-01	1.20e+04	5.48e-01	2.00e+00
3.00e+00	3.39e-01	1.25e-01	2.20e+04	3.49e-01	2.00e+00
4.00e+00	4.07e-01	6.25e-02	9.20e+04	3.06e-01	1.47e+00
5.00e+00	4.66e-01	3.12e-02	1.52e+05	3.09e-01	1.14e+00
6.00e+00	5.16e-01	1.56e-02	2.02e+05	3.09e-01	9.27e-01
7.00e+00	5.58e-01	7.81e-03	2.42e+05	3.13e-01	7.94e-01
8.00e+00	5.96e-01	3.91e-03	2.72e+05	3.10e-01	7.17e-01
9.00e+00	6.31e-01	1.95e-03	3.02e+05	3.02e-01	6.47e-01
1.00e+01	6.61e-01	9.77e-04	3.22e+05	3.01e-01	6.15e-01
1.10e+01	6.89e-01	4.88e-04	3.62e+05	3.12e-01	5.27e-01
1.20e+01	7.16e-01	2.44e-04	3.82e+05	3.05e-01	5.01e-01
1.30e+01	7.39e-01	1.22e-04	4.12e+05	3.28e-01	4.52e-01
1.40e+01	7.59e-01	6.10e-05	4.32e+05	3.16e-01	4.29e-01
1.50e+01	7.78e-01	3.05e-05	4.52e+05	3.05e-01	4.08e-01
1.60e+01	7.96e-01	1.53e-05	4.82e+05	3.11e-01	3.68e-01
1.70e+01	8.12e-01	7.63e-06	4.92e+05	3.01e-01	3.68e-01
1.80e+01	8.26e-01	3.81e-06	5.22e+05	3.09e-01	3.32e-01
1.90e+01	8.40e-01	1.91e-06	5.42e+05	3.08e-01	3.16e-01
2.00e+01	8.52e-01	9.54e-07	5.62e+05	3.08e-01	3.00e-01
2.10e+01	8.63e-01	4.77e-07	5.72e+05	3.02e-01	3.00e-01
2.20e+01	8.74e-01	2.38e-07	6.02e+05	3.20e-01	2.71e-01
2.30e+01	8.83e-01	1.19e-07	6.22e+05	3.12e-01	2.57e-01
2.40e+01	8.92e-01	5.96e-08	6.42e+05	3.17e-01	2.44e-01
2.50e+01	9.01e-01	2.98e-08	6.52e+05	3.00e-01	2.44e-01
2.60e+01	9.08e-01	1.49e-08	6.82e+05	3.15e-01	2.20e-01
2.70e+01	9.15e-01	7.45e-09	6.92e+05	3.05e-01	2.20e-01
2.80e+01	9.21e-01	3.73e-09	7.22e+05	3.18e-01	1.99e-01
2.90e+01	9.27e-01	1.86e-09	7.32e+05	3.05e-01	1.99e-01
3.00e+01	9.32e-01	9.31e-10	7.52e+05	3.07e-01	1.89e-01
3.10e+01	9.37e-01	4.66e-10	7.72e+05	3.21e-01	1.79e-01
3.20e+01	9.42e-01	2.33e-10	7.92e+05	3.19e-01	1.71e-01
3.30e+01	9.46e-01	1.16e-10	8.12e+05	3.16e-01	1.62e-01
3.40e+01	9.49e-01	5.82e-11	8.32e+05	3.33e-01	1.54e-01
3.50e+01	9.53e-01	2.91e-11	8.42e+05	3.19e-01	1.54e-01

Display

Print the true probability

fprintf('\nTrue probability \t %6.2e\n\n',Proba_true);
True probability 	 4.70e-11

Print the estimation and its confidence interval. This interval only in the asymptotic regime

fprintf('Lower bnd\t\tEstimate \t\tUpper bnd\n');
fprintf('%6.2e\t\t%6.2e \t\t%6.2e\n',Stat.interval(1),Proba_est,Stat.interval(2));
Lower bnd		Estimate 		Upper bnd
3.71e-11		5.10e-11 		6.32e-11

Show the last intermediate thresholds and estimated probability.

figure;
semilogy([Intern(end-10:end-1,2);threshold],[Intern(end-10:end-1,3);Proba_est],'-+');%,,'+');
hold on
semilogy(threshold*[1,1],[Stat.interval(1),Stat.interval(2)],'o-r');
semilogy(threshold,Proba_true,'^g')
legend('Estimate','Confidence','True')
grid on;
v = axis; v(2) =v(2)+0.01;
axis(v);
hold off
title('Test Estimator A')
xlabel('Threshold');
ylabel('Estimated');