%  test_ao_lpsd_variance_montecarlo
%
% Tests that the standard deviation returned by ao.dy in one
% frequency bin is equivalent to the matlab's std taking
% considering all realisations
%
% M Nofrarias 05-08-09
%
% $Id: test_ao_lpsd_variance_montecarlo.m,v 1.2 2009/08/11 14:20:10 miquel Exp $

% function test_ao_lpsd_variance_montecarlo()

clear

% data
nsecs = 500;
fs    = 5;
pl = plist('nsecs', nsecs, 'fs', fs, 'tsfcn', 'randn(size(t))');

%lpsd
pl2 = plist('Kdes', 20, 'Jdes', 100);

% loop
for  i = 1:100
  a(i) = ao(pl);
  b1(i) = lpsd(a(i),pl2);
end


%% compare

% select a realisation
index = 10;

% compare estimated std to the one using all realisations
clear rel
for i =1:len(b1(1))
  mn(i) = [mean(b1(:).y(i))];
  rel(:,i) = [std(b1(:).y(i))/mn(i) b1(index).dy(i)/b1(index).y(i)];
end

% k = find(b1(1).procinfo,'k')

figure
loglog(b1(1).x,rel')
hold
loglog(b1(1).x,abs(rel(2,:)-rel(1,:)),'-or')
ylabel('difference (%)')
% xlabel('Frequency [Hz]')
% ylabel('Rel. error')
% legend('montecarlo','estimated','difference')
% print -deps2c test_ao_lpsd_variance_montecarlo_1.eps

% end