%% SET the initial values for the Bayesian Anova.
load Data_INPUT
load PRIORS %priors for the parameters, set in Prior_Maker
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];


Pars_Selector=0; 
%set to 0 to set initial values of the scalar parameters to the modes of priors
%set to 1 for draws from (truncated) priors

%% SET THE INITIAL VALUES OF THE GRAND MEAN, ERROR VARIANCE, AND VARIANCES
% OF THE FACTOR LEVELS (IN THE REDUCED SPACE) TO THE MODES OF THE PRIORS:
if Pars_Selector==0
    INITIAL_VALS.GM=PRIORS.GM(1);
    INITIAL_VALS.sigma2=PRIORS.sigma2(2)/(PRIORS.sigma2(1)+1);
    
    
    for kk=1:1:N_factors
        %Initial Value for reduced space variances: modes of inverse gamma prior:
        tp=eval(['PRIORS.Red_space_var_',num2str(kk), '(2)/(PRIORS.Red_space_var_',num2str(kk),'(1)+1)']);
        INITIAL_VALS=setfield(INITIAL_VALS, ['Red_space_var_', num2str(kk)], tp);
    end
    
    
elseif Pars_Selector==1;
    
    %Initial Value for Grand Mean: draw from the normal prior.
    INITIAL_VALS.GM=PRIORS.GM(1)+sqrt(PRIORS.GM(2))*randn(1);
    
    %Initial Value for sigma2: draw from inverse gamma prior is likely a bad
    %ideas, due to the very large possible values. So truncate to below some
    %value:
    cutt=8;
    dummy=0;
    while dummy==0
        p1=1./randraw('gamma', [0,1/PRIORS.sigma2(2), PRIORS.sigma2(1)], 1);
        if p1<cutt
            INITIAL_VALS.sigma2=p1;
            dummy=1;
        end
    end
    
    
    for kk=1:1:N_factors
        %Initial Value for reduced space variances: : draw from inverse gamma prior truncated to
        %less than some cut off value:
        cutt=10;
        dummy=0;
        while dummy==0
            pars=eval(['PRIORS.Red_space_var_', num2str(kk)]);
            p1=1./randraw('gamma', [0,1/pars(2), pars(1)], 1);
            if p1<cutt
                INITIAL_VALS=setfield(INITIAL_VALS, ['Red_space_var_', num2str(kk)], p1);
                dummy=1;
            end
        end
    end
    
end

%% ALSO need initial values for the factor levels. 
%several choices:
Factor_selector=0;
% = 0: estimate levels using available data, set to these. 
% = 1: Draw in the reduced space using the initial values of the variances,
% hit with Q_M matrices. 

% Estimate the factor levels:
GM_est=nanmean(Data_INPUT(:));
f_levs_STACK=[];
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

if Factor_selector==0
    for kk=1:1:N_factors
        p1=f_levs_STACK(CSNL_Ext(kk)+1:1:CSNL_Ext(kk+1));
        INITIAL_VALS=setfield(INITIAL_VALS, ['Levels_', num2str(kk)], p1);
    end
    
elseif Factor_selector==0
    for kk=1:1:N_factors
        p1=eval(['PRIORS.Q_M_', num2str(kk)])*eval(['INITIAL_VALS.Red_space_var_', num2str(kk)])^(1/2)*randn(N_levels(kk)-1,1);
        INITIAL_VALS=setfield(INITIAL_VALS, ['Levels_', num2str(kk)], p1);
    end
    
end
    
%% FINALLY, need initial values for the missing observations in Data_INPUT:
% THESE ARE NOT NEEDED!
% will be drawn in the first iteration of the sampler, using the initial
% values of everything else.

%% SAVE IT
save INITIAL_VALS INITIAL_VALS