%% MAKE THE PRIORS FOR THE BAYESIAN ANOVA WITH MISSING VALUES
%NEED priors for:
% Grand Mean
% Error variance sig^2
% Variance for each factor. 
% NOTE THAT THE PRIORS FOR THE FACTORS ARE MEAN ZERO MVNs, with the
% specified covariance structure.

%% Read in the data matrix:
load Data_INPUT
%Data_INPUT=permute(Data_INPUT, [1,3,2]);
%get the number of factors, and the number of levels for each:
N_levels=size(Data_INPUT);
N_factors=length(N_levels);

%It will also be convenient to have the cumulative sum of the elements of
%N_levels:
CSNL=cumsum(N_levels);
CSNL_Ext=[0, CSNL];

%% Estimate the grand mean as the mean of all observed values. 
GM_est=nanmean(Data_INPUT(:));

%% estimate the factor levels:
% and the variance of the levels within each factor

f_levs_STACK=[];
vars_of_levels=NaN(N_factors,1);
for kk=1:1:N_factors
    dats=permute(Data_INPUT,[kk, setdiff([1:1:N_factors], kk)]);
    fac_est=nanmean(squeeze(nanmean(dats,2)),2)-GM_est;
    f_levs_STACK=[f_levs_STACK; fac_est(:)];
    vars_of_levels(kk)=var(fac_est(:));
end
%New factors begin at the entries of f_levs_STACK given by
%CSNL_Ext(1:1:N_levels-1)+1


%plot the truth and the new ones:
% load TRUE_FLS
% figure(1), clf
% for kk=1:1:N_factors
%     subplot(N_factors,1,kk)
%     plot(1:1:N_levels(kk), eval(['TRUE_FLS.Levels_', num2str(kk)]), 'ko-', 'linewidth', 2, 'markerfacecolor', 'k'), hold on
%     plot(1:1:N_levels(kk), f_levs_STACK(CSNL_Ext(kk)+1:1:CSNL_Ext(kk+1)), 'r+--', 'linewidth', 1.5), hold on
% end

%% First construct arrays that correspond to the levels of each factor:
FAC_ARRAY=[];

for jj=1:1:N_factors
    fac_arr=permute(NaN(N_levels), [jj, setdiff([1:1:N_factors], jj)]);
    %fill it
    for kk=1:1:N_levels(jj)
        fac_arr(kk,:,:)=f_levs_STACK(kk+sum(N_levels)-sum(N_levels(jj:1:end)));
    end
    %undo the permutation:
    [vals, ord]=sort([jj, setdiff([1:1:N_factors], jj)]);
    fac_arr=permute(fac_arr, ord);
    FAC_ARRAY=cat(N_factors+1,FAC_ARRAY,fac_arr);
end


%% Estimate the residual between each observation and the corresponding
% combination of the factors. 
differs=Data_INPUT - (squeeze(sum(FAC_ARRAY,N_factors+1))+GM_est);
Est_Var=nanvar(differs(:));

%% Prior for the grandmean:
%Mean: average of all available observations.
%Var: 16*variance of residuals, which are estimated by remomving the
%estimated grand mean and factor levels from the observations. 
%IDEA: very vague (high var), but centered at a reasonable point.
PRIORS.GM=[GM_est, 16*Est_Var]; 

%% Prior for the variance of the errors
% Set it to be the equivalent of one observation, with average squared
% deviation equal to the variance of the estimated residuals. 
PRIORS.sigma2=[0.5,Est_Var/2];
% Note, caliing the pars a and b, and the pars of the scaled inverse chi^2
% nu and s^2, we have nu=2a and s^2=b/a. These give, in the bayesian sense,
% the number of prior samples and the mean squared deviation.

%% PRIORS for the reduced space variances (one for each factor):
%Set each to correspond to half an observation (RECALL there might not be
%many levels for some factors, so we do not want to unduly influence the
%inference) with average squared deviation given by the variance of the
%estimated levels.
for jj=1:1:N_factors
    p_pars=[0.25, vars_of_levels(jj)/4];
    PRIORS=setfield(PRIORS, ['Red_space_var_', num2str(jj)], p_pars);
end

%% WE will also need the transform matrices, Q_M in the notes, for each of
% the factors

for jj=1:1:N_factors
    PRIORS=setfield(PRIORS, ['Q_M_', num2str(jj)], T_mat_Genner(N_levels(jj)));
end


%% 

save PRIORS PRIORS
