% quantile averaging in the case where some subjects have few errors
% if in any condition any subject has fewer than 11 errors, this error RT
% distribution is ignored.
% NB The lines below are exclusively for quantile fitting!

dmat_setversion(.31)

min_RT = 0.3; 
max_RT = 2.5; 
q      = [.1 .3 .5 .7 .9]; % the quantiles
N_cond = 8;
N_sub  = 2; %number of subjects; for this illustration we'll just duplicate
            %the data from "pp14.mat" twice  
ini    = cell(1,N_cond);
data   = struct('CorrectObs', ini,   ...
                'IncorrectObs', ini, ...
                'CorrectEdges', ini, ...
                'IncorrectEdges', ini);
            
for i_sub = 1:N_sub
    if i_sub==1
        fn = 'C:\Documents and Settings\Jepma\Bureaublad\difussiedata\pp14.mat';
    end
    if i_sub==2
        fn = 'C:\Documents and Settings\Jepma\Bureaublad\difussiedata\pp14.mat';
    end
    
    s  = load(fn); % All data of one Ss are now in s.data 
    %(1st column = condition, 2nd column = accuracy, 3rd column = RT)
    % now we need to calculate stuff for each condition seperately.

    for i = 1:N_cond        
        cond_index = find(s.data(:,1)==i & s.data(:,3)>min_RT & s.data(:,3)<max_RT); 
         % index contains row numbers corresponding to conditions
        
        correct_index     = intersect(cond_index, find(s.data(:,2)==1));
        total_correct     = length(correct_index);

        incorrect_index   = intersect(cond_index, find(s.data(:,2)==0));
        total_incorrect   = length(incorrect_index);
        
        % determine whether each condition has > 10 observations:
        enough_corrects(i_sub, i)   = total_correct > 10;        
        enough_incorrects(i_sub, i) = total_incorrect > 10;        
    end    
end
% Now we've determined which conditions have "enough" observations to be
% taken into account.
sum_cor   = sum(enough_corrects,1);
sum_incor = sum(enough_incorrects,1);

% Actual averaging:
for i_sub = 1:N_sub
    if i_sub==1
        fn = 'C:\Documents and Settings\Jepma\Bureaublad\difussiedata\pp14.mat';
    end
    if i_sub==2
        fn = 'C:\Documents and Settings\Jepma\Bureaublad\difussiedata\pp14.mat';
    end
    
    s  = load(fn); % All data of one Ss are now in s.data

    for i = 1:N_cond        
        cond_index        = find(s.data(:,1)==i & s.data(:,3)>min_RT & s.data(:,3)<max_RT); 
        
        correct_index     = intersect(cond_index, find(s.data(:,2)==1));
        total_correct     = length(correct_index);

        incorrect_index   = intersect(cond_index, find(s.data(:,2)==0));
        total_incorrect   = length(incorrect_index);
                
        if (enough_corrects(i_sub,i) == 1)
            if ( length(data(i).CorrectEdges) == 0 ) %no data yet
                data(i).CorrectEdges   = (quantile(s.data(correct_index,3), q)) / sum_cor(i); 
                data(i).CorrectObs     = diff([0 q 1]) * total_correct;
            else %data already exist
                data(i).CorrectEdges   = data(i).CorrectEdges + (quantile(s.data(correct_index,3), q)) / sum_cor(i); 
                data(i).CorrectObs     = data(i).CorrectObs + diff([0 q 1]) * total_correct;
            end                
        end

        if (enough_incorrects(i_sub,i) == 1)            
            if ( length(data(i).IncorrectEdges) == 0 ) %no data yet
                data(i).IncorrectEdges = (quantile(s.data(incorrect_index,3), q)) / sum_incor(i);
                data(i).IncorrectObs   = diff([0 q 1]) * total_incorrect;
            else %data already exist
                data(i).IncorrectEdges = data(i).IncorrectEdges + (quantile(s.data(incorrect_index,3), q)) / sum_incor(i);                
                data(i).IncorrectObs   = data(i).IncorrectObs + diff([0 q 1]) * total_incorrect;            
            end
        end        
    end    
end


% Now back to regular fitting commands:
%% this example fits the All free model
% The designmatrix for the Ter-model is {[1,0;1,0;1,0;1,0;0,1;0,1;0,1;0,1],[1,0;1,0;0,1;0,1;1,0;1,0;0,1;0,1],[1,0;0,1;1,0;0,1;1,0;0,1;1,0;0,1],[1;1;1;1;1;1;1;1],[1,0;1,0;1,0;1,0;0,1;0,1;0,1;0,1],[1;1;1;1;1;1;1;1],[1,0;0,1;1,0;0,1;1,0;0,1;1,0;0,1]}
% The designmatrix for the a-model is {[1,0,0,0;1,0,0,0;0,1,0,0;0,1,0,0;0,0,1,0;0,0,1,0;0,0,0,1;0,0,0,1],[1;1;1;1;1;1;1;1],[1,0;0,1;1,0;0,1;1,0;0,1;1,0;0,1],[1;1;1;1;1;1;1;1],[1,0,0,0;1,0,0,0;0,1,0,0;0,1,0,0;0,0,1,0;0,0,1,0;0,0,0,1;0,0,0,1],[1;1;1;1;1;1;1;1],[1,0;0,1;1,0;0,1;1,0;0,1;1,0;0,1]}
% The designmatrix for the v-model is {[1,0;1,0;1,0;1,0;0,1;0,1;0,1;0,1],[1;1;1;1;1;1;1;1],[1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1;1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1],[1;1;1;1;1;1;1;1],[1,0;1,0;1,0;1,0;0,1;0,1;0,1;0,1],[1;1;1;1;1;1;1;1],[1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1;1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1]}

options(1) = struct('DesignMatrix',{{[1,0,0,0;1,0,0,0;0,1,0,0;0,1,0,0;0,0,1,0;0,0,1,0;0,0,0,1;0,0,0,1],[1,0;1,0;0,1;0,1;1,0;1,0;0,1;0,1],[1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1;1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1],[1;1;1;1;1;1;1;1],[1,0,0,0;1,0,0,0;0,1,0,0;0,1,0,0;0,0,1,0;0,0,1,0;0,0,0,1;0,0,0,1],[1;1;1;1;1;1;1;1],[1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1;1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1]}},...
   'Display','off',...
   'EWMA',struct('l',0.01,...
   'L',1.5,...
   's',0.5),...
   'EstimationMethodScalar',6,...
   'FixedBinEdges',[0.3,0.36,0.42,0.52,0.8;0.38,0.47,0.56,0.7,1],...
   'FixedValues',[],...
   'Guess',[],...
   'GuessMethodScalar',1,...
   'LongSimplexRuns',1,...
   'MaxIter',5000,...
   'Name','all_free',...
   'NoFitting',0,...
   'NonparametricBootstrap',0,...
   'ObjectiveDecimals',7,...
   'OutlierMax',2.5,...
   'OutlierMin',0.3,...
   'OutlierTreatment','Absolute cut-off',...
   'ParameterDecimals',7,...
   'ParametricBootstrap',0,...
   'Percentiles',[10,30,50,70,90;10,30,50,70,90],...
   'ShortSimplexRuns',3,...
   'SpecificBias',[]);
%% Fit model
output = runqueue_q(data,options);  %NB runqueue_q fits quantiles!

for a = 1:8,
        [cdf pxis1(a)] = cdfdif(0,0,output.Minimum(a,:));
end

%% Show a small summary table
qtable(output);
disp('      Queue finished. Check output structure array for results.');

%% Print output to file
fid  = fopen('group_level_AllFreeModel.out', 'w'); % create a file and open it for writing. 

fprintf(fid, 'The name of the data file is:\n');
fprintf(fid, '%s\n', fn);
fprintf(fid, 'The number of free parameters is:\n');
fprintf(fid, '%1i\n', output.Df);
fprintf(fid, 'The BIC value is:\n');
fprintf(fid, '%3.4f\n', output.FitInfo.BIC);
fprintf(fid, 'The best parameter estimates are:\n');
fprintf(fid, '%3.4f\n', output.Minimum);
fprintf(fid, 'The standard errors are:\n');
fprintf(fid, '%3.4f\n', output.StdErr);

fclose(fid); %closes the file.