%% MAKE SOME TEST DATA TO TEST BAYESIAN ANOVA
% Idea: make code general, so it will work with an K dimensional input
% array. 
% Assume: Any straight line through the data will hit at least one observed
% value. (In other words, if all but one of the factor levels are held
% constant, there is at an observation at at least one of the level of the
% last factor.)

% For this example, K=2

%% Set the three vectors of factors:
%number of factors:
NN=2;
%number of levels:
N1=20;
N2=15;

% set the std for each, IN THE REDUCED SPACE
sig_red1=10;
sig_red2=15;

% get the transform matrices for each factor (different for each on account
% of the different number of levels)
T_Mat1=T_mat_Genner(N1);
T_Mat2=T_mat_Genner(N2);

% to get the factors in the full space, generate in the reduced, and hit
% with the transfrom matrices:

fac1=T_Mat1*sig_red1*randn(N1-1,1);
fac2=T_Mat2*sig_red2*randn(N2-1,1);
sum(fac1)
sum(fac2)

%Set the grand mean and the standard deviation of the errors:
GM=5;
Var_errs=10^2;

%% PLOT the factors:
figure(1), clf
subplot(NN,1,1)
plot(1:1:N1, fac1, 'ko-', 'linewidth', 2, 'markerfacecolor', 'k'), hold on
subplot(NN,1,2)
plot(1:1:N2, fac2, 'ko-', 'linewidth', 2, 'markerfacecolor', 'k'), hold on

%% MAKE the data array:
%first make an array corresponding to each factor:
fac1_array=ones(N1, N2);
for kk=1:1:N1
    fac1_array(kk,:)=fac1(kk);
end

fac2_array=ones(N1, N2);
for kk=1:1:N2
    fac2_array(:,kk)=fac2(kk);
end


Data_Array=fac1_array+fac2_array+GM+sqrt(Var_errs)*randn(N1,N2);

%for ii=1:1:N1
%     for jj=1:1:N2
%         for kk=1:1:N3
%            Data_Array(ii,jj,kk)=fac1(ii)+fac2(jj)+fac3(kk)+GM+sqrt(Var_errs)*randn(1,1);
%         end
%     end
% end
%% Estimate the grand mean, and the factor levels:
GM_est=mean(Data_Array(:))

fac1_est=mean(Data_Array,2)-GM_est;
fac1_est=fac1_est(:);

fac2_est=mean(Data_Array,1)-GM_est;
fac2_est=fac2_est(:);

%plot the estimated factor levels:
figure(1), gcf
subplot(NN,1,1)
plot(1:1:N1, fac1_est, 'r+--', 'linewidth', 1.5), hold on
subplot(NN,1,2)
plot(1:1:N2, fac2_est, 'r+--', 'linewidth', 1.5), hold on

%% NOW knock out half of the values, re-estimate the grand mean and factor
% levels:
% DO IT by knocking out all entries corresponding to index n-tuples with
% odd sums:
Data_Array_N=Data_Array;
for ii=1:1:N1
    for jj=1:1:N2
           if floor(sum(ii+jj)/2)~=sum(ii+jj)/2
               Data_Array_N(ii,jj)=NaN(1);
           end
    end
end


%% RE-ESTIMATE the grand mean, and factor levels:
GM_est_N=nanmean(Data_Array(:))

fac1_est_N=nanmean(Data_Array_N,2)-GM_est_N;
fac1_est_N=fac1_est_N(:);

fac2_est_N=nanmean(Data_Array_N,1)-GM_est_N;
fac2_est_N=fac2_est_N(:);


%plot the estimated factor levels:
figure(1), gcf
subplot(NN,1,1)
plot(1:1:N1, fac1_est_N, 'b+--', 'linewidth', 1.5), hold on
subplot(NN,1,2)
plot(1:1:N2, fac2_est_N, 'b+--', 'linewidth', 1.5), hold on

subplot(NN,1,1)
title('BLACK: True. Red: Estimate, all data. Blue: Estimate, half NaNs')


%% SAVE a bunch of stuff:
%the complete data array:
Data_Array_Complete=Data_Array;
save Data_Array_Complete Data_Array_Complete

%With half repalced by NaNs
Data_Array_wNANs=Data_Array_N;
save Data_Array_wNANs Data_Array_wNANs
Data_INPUT=Data_Array_wNANs;
save Data_INPUT Data_INPUT

%A VECTOR composed of the FACTORS, in order:
% Factor_Vector=[fac1(:); fac2(:); fac3(:)];
% save Factor_Vector Factor_Vector
% ACTUALLY prefer if this matches up with how the factor levels are treated
% in the main code:
TRUE_FLS=struct('Levels_1', []);
for kk=1:1:NN
    TRUE_FLS=setfield(TRUE_FLS, ['Levels_', num2str(kk)], eval(['fac',num2str(kk)]));
end
save TRUE_FLS TRUE_FLS


%A VECTOR composed of the (reduced space) standard deviations, in
%order:
% Factors_ReducedSpace_Stds=[sig_red1; sig_red2; sig_red3];
% save Factors_ReducedSpace_Stds Factors_ReducedSpace_Stds
% %noise error:
% Noise_Var=Var_errs;
% save Noise_Var Noise_Var
% %grand mean:
% Grand_Mean=GM;
% save Grand_Mean Grand_Mean
% ACTUALLY prefer if this matches up with how the factor levels are treated
% in the main code:
TRUE_PARS.GM=GM;
TRUE_PARS.sigma2=Var_errs;
for kk=1:1:NN
    TRUE_PARS=setfield(TRUE_PARS, ['Red_space_var_', num2str(kk)], eval(['sig_red',num2str(kk), '^2']));
end

save TRUE_PARS TRUE_PARS