Test for estimator_B
Estimation of a weak probability
Authors
- Implementation: Teddy Furon
- Science: Arnaud Guyader, Nicolas Hengartner, Eric Matzner-Lober
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
- x a white gaussian noise in R^20
- lies in the hypercone of axis (1,0,...,0) and angle acos(0.95)
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
