%% Confirmatory sampling simulation to mirror 12 datasets.
% Specify inputs

clear all

% means [S, T, H]
mu = [15 15 15];

% Standard deviations [S, T, H]
s = [ 3 ; 3 ; 3 ];

%n of samples
n = 100;

%n of simulations (datasets)
nSims = 100;

% Input Corrs 
inputCorr_ST = .5;
inputCorr_SH = .5;
inputCorr_TH = .5;
    
  %correlation matrix [S, T, H]
    R = [ 1.0,          inputCorr_ST, inputCorr_SH,
          inputCorr_ST, 1.0,          inputCorr_TH,
          inputCorr_SH, inputCorr_TH, 1.0 ];
      
%% generate samples
      
V = s*s';
r = {};

for i = 1:nSims
% Co-variance matrix:
covariance = V.*R;

% Generate Samples
r{i,1} = mvnrnd( mu, covariance, n ); %cell containing all samples
r{i,1} = round(r{i,1}); %round to whole numbers

% View descriptives
maxVal{i} = max(r{i});
minVal{i} = min(r{i});
medians{i} = median(r{i});
means{i} = mean(r{i});
stdevs{i} = std(r{i});
corrs{i} = corrcoef(r{i});

end %for multivariate

%% do analyses within each of 100 sims (within dataset analysis)

%define cells for storage over iteration
S = {};
T = {};
H = {};
model_TonS = {};
predicted_T = {};
Rt = {};
model_SonT = {};
predicted_S = {};
Rs = {};

%loop over n (100) datasets
for i = 1:nSims %for i

S{i,1} = r{i,1}(:,1);
T{i,1} = r{i,1}(:,2);
H{i,1} = r{i,1}(:,3);
    
model_TonS{i,1} = polyfit(S{i},T{i},1); %get model fit predicting T from S
predicted_T{i,1} = polyval(model_TonS{i,1},S{i}); %use model to predict T

Rt{i} = (T{i,1} - predicted_T{i,1}) * -1; %get residuals Rt %must reverse score T

model_SonT{i,1} = polyfit(T{i},S{i},1); %get model fit predicting S from T
predicted_S{i,1} = polyval(model_SonT{i,1},T{i}); %use model to predict S

Rs{i} = S{i,1} - predicted_S{i,1}; %get residuals Rs

corr_Rs_S(i,1) = corr(Rs{i},S{i}); %criterion corr: r(S,Rs)
corr_Rs_T(i,1) = corr(Rs{i},T{i}); %should always be 0 (independent)

corr_S_T(i,1) = corr(S{i},T{i}); %accuracy corr

two_resids_corr(i) = corr(Rt{i},Rs{i}); %correlate both residuals OVER 100 datasets

end %for i

acc_resid_corr = corr(corr_Rs_S,corr_S_T); %confirmation value in paper: r(r(Rs,S),r(S,T)

%% stack cells to do analyses over all nSims (over all 100 dataset analysis)

stackedSims = cell2mat(r); %stack all 100 datasets on top of each other

S = stackedSims(:,1); %turn into vector of S
T = stackedSims(:,2); %vector of T
H = stackedSims(:,3); %vector of H

model_TonS = polyfit(S,T,1); %get model fit predicting T from S
predicted_T = polyval(model_TonS,S); %use model to predict T

Rt = (T - predicted_T) * -1; %get residuals Rt %must reverse score

model_SonT = polyfit(T,S,1); %get model fit predicting S from T
predicted_S = polyval(model_SonT,T); %use model to predict S

Rs = S - predicted_S; %get residuals Rs

two_resids_corr = corr(Rt,Rs); %correlate both residuals over 100000 observations

S_T = (S - T); %compute S-T

abs_S_T = abs(S_T); %compute |S-T|
abs_Rs = abs(Rs); %compute |Rs|

corr(S,Rs); %target correlation r(Rs,S)
corr(T,Rs); %should always be 0 (independent)

