% Particle Filter Demo MAE
% 7/6/20 -- Method 7 uses 50% whiskers, 1 equilibration
% 1/16/19 -- Replaces 100% whiskers with 25% whiskers
% 1/4/19 -- Tabulates MAE
% New Particle Filter Demo 2/16/18
% Double Whiskers 3/9/18
% Beta scatter added 5/12/18
% Single Whisker added 10/13/18
% Particulate Anticipation added 10/30/18
% y(t) = x(t) + epsilon(t)
% x(t) = x(t-1) + eta(t)
% N = # of particles.  Must be even for Mixed Particulate Adaptive
% equilibration
clear
bugtrap = 1
Reps = 1000  % Reps <= 1000
T = 100
N = 100   % N = 100, 1000, 10000, 100000
alpha = 1.7
ceps = 1
ceta = .25
betaeps = .3
betaeta = 0
param = 1 
%  = 0 for continuous zeta parameterization
%  = 1 for mean-focus mu parameterization

% adaptive 
% = 0 for basic
% = 1 for fully adaptive
% = 2 for bounded
% = 3 for particulate with fraction thetaAdaptive adaptive.  
% thetaAdaptive = .5 for half Adaptive Whiskers;
% equilibrate 
% = 0 for no equilibration with fractional particulate adaptive
% = 1 for single equiibration
% = 2 for double equilibration
whiskers = 0
% = 0 for basic
% = 1 for single whisker (new)
% = 2 for multiple whiskers (old)
scatter = 0
% = 0 for basic
% = 1 for Beta(a,a) scatter
continuous = 1 
% = 0 for step draw resampling
% = 1 for continuous draw resampling
% random 
% = 0 for stratified sampling, resampling in basic
% = 1 for random 

seed = 1
rng(seed)
% Generate and plot random data
% Uses Nolan RNG etc.  Defaults to 0 parameterization (continuous)
tic
eps = stablernd(1000*T, [alpha betaeps ceps 0], param);
eps = reshape(eps, 1000, T);
eps = eps(1:Reps, :);
'eps done'
eta = stablernd(1000*T, [alpha betaeta ceta 0], param);
eta = reshape(eta, 1000, T);
eta = eta(1:Reps, :);
'eta done'
x = cumsum(eta')';
y = x + eps;
'x y done'
if bugtrap
  MAeps = mean(abs(eps(:,11:T)));
  avgMAeps = mean(MAeps)
  medMAeps = median(MAeps)
  Q90MAeps = quantile(MAeps, .9)
  maxMAeps = max(MAeps)
end % if bugtrap
toc
% ~ Reps x T

% 1: Basic Random
% 2: Basic Stratified
% 3: Adaptive
% 4: 50% Whisker-Adaptive
% 5: 25% Basic/Whisker-Adaptive
% 6: 25%, 1 equilibration
% 7: 50%, 1 equilibration

tic
ESSbar = nan(7, Reps);
ESSmin = nan(7, Reps);
MSE = nan(7, Reps);

% j = 1
% %  Basic with random sampling
% adaptive = 0;
% random = 1;
% equilibrate = 0;
% for r = 1:Reps
%   r = r  
% [ xF, pF, logl, ESSbar(j,r), ESSmin(j,r), MAE(j,r)] =  ... 
%    ParticleFilterDemoFunction(y(r,:), x(r,:), alpha, ceps, ceta, betaeps, N, ... 
%      param, adaptive, equilibrate, random);
% end % for r
% 'basic random done'
% toc
% 
% j = 2
% %  Basic with stratified sampling
% random = 0
% for r = 1:Reps
%   r = r  
% [xF, pF, logl, ESSbar(j,r), ESSmin(j,r), MAE(j,r)] =  ... 
%    ParticleFilterDemoFunction(y(r,:), x(r,:), alpha, ceps, ceta, betaeps, N, ... 
%      param, adaptive, equilibrate, random);
% end % for 
% 
% j = 3
% % Simple Adaptive
% adaptive = 1
% for r = 1:Reps
%   r = r  
% [xF, pF, logl, ESSbar(j,r), ESSmin(j,r), MAE(j,r)] =  ... 
%    ParticleFilterDemoFunction(y(r,:), x(r,:), alpha, ceps, ceta, betaeps, N, ... 
%      param, adaptive, equilibrate, random);
% end % for 
% 
% j = 4
% % 50% Whisker-Adaptive
% adaptive = 3
% thetaAdaptive = .50
% equilibrate = 0
% random = 0
% for r = 1:Reps
%   r = r  
% [xF, pF, logl, ESSbar(j,r), ESSmin(j,r), MAE(j,r)] =  ... 
%    ParticleFilterDemoFunction(y(r,:), x(r,:), alpha, ceps, ceta, betaeps, N, ... 
%      param, adaptive, equilibrate, random, thetaAdaptive);
% end % for r 
% 
% j = 5;
% % 25% Whisker-Adaptive
% adaptive = 3
% thetaAdaptive = .25
% equilibrate = 0
% for r = 1:Reps
%   r = r  
% [xF, pF, logl, ESSbar(j,r), ESSmin(j,r), MAE(j,r)] =  ... 
%    ParticleFilterDemoFunction(y(r,:), x(r,:), alpha, ceps, ceta, betaeps, N, ... 
%      param, adaptive, equilibrate, random, thetaAdaptive);
% end % for r 
% 
% j = 6;
% % Basic + W-A, 1 equilibration
% adaptive = 3
% thetaAdaptive = .25
% equilibrate = 1
% for r = 1:Reps
%   r = r  
% [xF, pF, logl, ESSbar(j,r), ESSmin(j,r), MAE(j,r)] =  ... 
%    ParticleFilterDemoFunction(y(r,:), x(r,:), alpha, ceps, ceta, betaeps, N, ... 
%      param, adaptive, equilibrate, random, thetaAdaptive);
% end % for r 

j = 7;
% Basic + W-A, 1 equilibration
random = 0
adaptive = 3
equilibrate = 1
thetaAdaptive = .5
for r = 1:Reps
  r = r  
[xF, pF, logl, ESSbar(j,r), ESSmin(j,r), MAE(j,r)] =  ... 
   ParticleFilterDemoFunction(y(r,:), x(r,:), alpha, ceps, ceta, betaeps, N, ... 
     param, adaptive, equilibrate, random, thetaAdaptive);
end % for r
'Reps done'

toc

figure
plot(1:Reps, ESSbar(1:j,:), 'linewidth', 2)
title('ESSbar')
legend('Basic Random', 'Basic Stratified', 'Adaptive', '50% Whisker-Adaptive', ...
  '25% Basic/Whisker-Adaptive', ...
  'B/W-A, 1 equilibration', 'B/W-A, 2 equilibrations')
xlabel('replication')
grid on
figure
plot(1:Reps, ESSmin(1:j,:), 'linewidth', 2)
title('ESSmin')
legend('Basic Random', 'Basic Stratified', 'Adaptive', '50% Whisker-Adaptive', ...
  '25% Basic/Whisker-Adaptive', ...
  'B/W-A, 1 equilibration', 'B/W-A, 2 equilibrations')
xlabel('replication')
grid on
figure
plot(1:Reps, MAE(1:j,:), 'linewidth', 2)
title('MAE')
legend('Basic Random', 'Basic Stratified', 'Adaptive', '50%Whisker-Adaptive', ...
  '25% Basic/Whisker-Adaptive', ...
  'B/W-A, 1 equilibration', 'B/W-A, 2 equilibrations')
xlabel('replication')
grid on
'plots done'
toc
alpha = alpha
N=N
Reps=Reps
(1:7) + .00000001
avgMAE = mean(MAE(1:j,:)')
medMAE = median(MAE(1:j,:)')
Q90MAE = quantile(MAE', .90)
maxMAE = max(MAE(1:j,:)')
avgESSmin = mean(ESSmin(1:j,:)')
minESSmin = min(ESSmin(1:j,:)')
Q10ESSmin = quantile(ESSmin',.10)
medESSmin = median(ESSmin(1:j,:)')
'1: Basic Random'
'2: Basic Stratified  '
'3: Adaptive              '
'4: 50% Whisker-Adaptive  '
'5: 25% Basic/Whisker-Adaptive  '
'6: 1 equilibration  '
'7: 2 equilibrations ' 

return


