function calculate_error_dist_given_report(model_index)
% this function loads model fits for each patient and
% outputs patient_pain_reports
load('patient_data_all.mat');
% load modeling frameworks
load('model_names.mat');
modelname = model_names{model_index};
num_patients = length(patient_IDs);
std_res = zeros(num_patients,1);

for ii=1:num_patients
    % load patient data into convenient form (just useful info on specified time)
    filename = patient_IDs{ii}
    
    % read in model report for each patient to assess degree of freedom
    report = csvread(strcat('Patient_fixk0_',modelname,'/',filename,'_report_fixk0.csv'));

    % plot distribution of errors (residuals)
    residuals = report(2,:)-report(3,:);
%     figure
%     hist(residuals)
    % print std of residuals
    std_res(ii) = std(residuals);
    
    % test residuals for normality
%     x = (residuals-mean(residuals))/std(residuals);
%     [h,p] = kstest(x);
%     p
end

mean(std_res)
median(std_res)

end


% TAKE AWAY: 2 OF THE 39 PATIENTS HAVE SIGNIFICANTLY NON-NORMAL ERRORS