Test for estimator_B

Estimation of a weak probability

Authors

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

Contents

Example

What is the probability that

Data

dim = 20;
threshold = 0.95;
n = 200; % Total number of particles
t = 30;
% 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_B(n,'P',Proba_true)
Stat_pre = 

      NbIter: 4745
    interval: [2x1 double]
    rel_bias: 0
     rel_var: 0.1263

Estimation

[Proba_est, Stat, Intern] = estimator_B(threshold,dim,'n',n,'N_sample',200,'t',t);
Lower = Stat.interval(1);
Upper = Stat.interval(2);
List of all arguments of estimator_B
            a: 0.9500
          dim: 20
     GENERATE: @GENERATE_INTERNAL
     max_iter: 50000
      minheap: 0
       MODIFY: @MODIFY_INTERNAL
           mu: 2
     mu_decay: 0.0500
            n: 200
     N_sample: 200
         rate: 0.1000
        SCORE: @SCORE_INTERNAL
            t: 30
    threshold: 0.9500
      verbose: 1

Iteration 	Threshold 	Prob 		mu 		Trials 
1		2.57e-04	1.00e+00	2.00e+00	2.00e+02
200		2.09e-01	3.69e-01	2.00e+00	6.17e+03
400		3.41e-01	1.35e-01	2.00e+00	1.22e+04
600		4.38e-01	4.97e-02	2.00e+00	1.82e+04
800		5.13e-01	1.82e-02	1.55e+00	2.42e+04
1000		5.79e-01	6.69e-03	1.26e+00	3.02e+04
1200		6.25e-01	2.45e-03	1.14e+00	3.62e+04
1400		6.71e-01	9.01e-04	8.80e-01	4.22e+04
1600		7.07e-01	3.30e-04	7.94e-01	4.82e+04
1800		7.43e-01	1.21e-04	6.81e-01	5.42e+04
2000		7.71e-01	4.45e-05	6.15e-01	6.02e+04
2200		7.95e-01	1.63e-05	5.84e-01	6.62e+04
2400		8.17e-01	5.99e-06	5.55e-01	7.22e+04
2600		8.35e-01	2.20e-06	4.76e-01	7.82e+04
2800		8.49e-01	8.07e-07	4.76e-01	8.42e+04
3000		8.66e-01	2.96e-07	4.29e-01	9.02e+04
3200		8.81e-01	1.09e-07	4.08e-01	9.62e+04
3400		8.92e-01	3.99e-08	3.87e-01	1.02e+05
3600		9.04e-01	1.46e-08	3.50e-01	1.08e+05
3800		9.14e-01	5.37e-09	3.32e-01	1.14e+05
4000		9.23e-01	1.97e-09	3.16e-01	1.20e+05
4200		9.31e-01	7.23e-10	2.85e-01	1.26e+05
4400		9.38e-01	2.65e-10	2.57e-01	1.32e+05
4600		9.44e-01	9.74e-11	2.32e-01	1.38e+05
4800		9.49e-01	3.57e-11	2.32e-01	1.44e+05

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',Lower,Proba_est,Upper);
Lower bnd		Estimate		Upper bnd
1.58e-11		3.15e-11		6.17e-11

Show the last intermediate thresholds and estimated probability.

figure;
subplot(1,2,1)
semilogy([Intern(end-10:end-1,2);threshold],[Intern(end-10:end-1,3);Proba_est],'-+');%,,'+');
hold on
semilogy(threshold*[1,1],[Lower,Upper],'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 B')
xlabel('Threshold');
ylabel('Estimated');

subplot(1,2,2)
x = linspace(1/2*Stat.interval(1),2*Stat.interval(2),100);
mu = -n*log(Proba_est)*log(1-1/n);
sig2 = -n*log(Proba_est)*log(1-1/n)^2;
pdf_estimator = @(x,m,s2) exp(-(log(x)-m).^2/2/s2)./x/sqrt(2*pi*s2);
%PDF of the estimator - log normal
plot(x,pdf_estimator(x,mu,sig2)); hold on;
x = [Proba_true,Stat.interval(1),Stat.interval(2),Proba_est];
y = pdf_estimator(x,mu,sig2);
stem(x(1),y(1),'^g')
stem(x(2:3),y(2:3),'r');
stem(x(4),y(4),'+')
hold off
xlabel('probability')
ylabel('pdf')
title('pdf of the estimator')
hold off