clear all; close all;

addpath ../../Toolboxes/Palamedes

%% Define what to analyse

% Which subjects?
subjects = [1:34];

% Which data?
whichdata = 'clean';  % alternative: 'raw'

if whichdata == 'clean'
    load('P:/3018029.09/data_sd_clean.mat');
    folder_str = 'Clean data compared to simulated data';
elseif whichdata == 'raw'
    load('P:/3018029.09/data_sd.mat');
    folder_str = 'Raw data compared to simulated data';
end
alldata = data_sd; clear data_sd

% Append stimulus difficulty to data
difficulty_levels = [6 5 4 3 2 1 1 2 3 4 5 6];

prevDifficulty = difficulty_levels(alldata.prevStimLevel)';
currDifficulty = difficulty_levels(alldata.currStimLevel)';

alldata = [alldata, table(prevDifficulty), table(currDifficulty)];

% Which simulated data?
load('P:/3018029.09/simulatedResponses_1000_iterations.mat')
iterations = 1000;

%Parameter grid defining parameter space through which to perform a
%brute-force search for values to be used as initial guesses in iterative
%parameter search.
searchGrid.alpha = -0.2:.001:0.2;
searchGrid.beta = logspace(0,3,101);
searchGrid.gamma = [0:.01:.06];
searchGrid.lambda = [0:.01:.06];

%% P-repeat: real data 

% RT

for rt = 1:2

for sj = unique(subjects)
    
    % individual sj data
    data = alldata(alldata.subject == subjects(sj), :);
        
    nbins = 6;
      
    for u = 1:nbins
                
        % first select trials for this bin
        trls = find(data.prevDifficulty == u);
                
        % find median rt
        datartsorted    = sortrows(data.prevRt(trls));
        medianrt        = datartsorted(floor(length(datartsorted)/2));
        
      % select low or high rts
        if rt == 1 
            trlssplit = trls(data.prevRt(trls) < medianrt);
        elseif rt == 2
            trlssplit = trls(data.prevRt(trls) >= medianrt);
        end
                
        % =========================================== %
        % history-dependent logistic
        % =========================================== %
        
        % previous response
        resps = [0 1];
        for r = 1:2
            
            % take a further subset of trials with previous response 1 or 0
            thisdat = data(trlssplit,:);
            thisdat = thisdat(thisdat.prevResponseDir == resps(r),:);
            
            [bias, slope, lapseLow, lapseHigh] = fitLogistic(...
                thisdat.currRelativeCoherence, thisdat.currResponseDir);
            
            % save betas
            grandavg{rt}.logisticHistory(sj, r, u, :) = [bias slope];
        end    
           
    end
          
end

% ========================================================= %
% transform history-dependent bias into repetition bias
% ========================================================= %

resp1 = -1 * squeeze(grandavg{rt}.logisticHistory(:, 1, :, 1)); % take the 1st beta weigth only, ignore history-dependent slope
resp2 =      squeeze(grandavg{rt}.logisticHistory(:, 2, :, 1));

% since this is centred at 0.5, treat it that way
grandavg{rt}.repetition = (resp1 + resp2) ./ 2;

% ========================================================= %
% transform from log(odds) to probability
% ========================================================= %

logOdds2Probability = @(x) exp(x) ./ (1 + exp(x));
grandavg{rt}.repetition = logOdds2Probability(grandavg{rt}.repetition);

end

% save means and SEs for prev low and high rt

mean_repetition_lowrt = mean(grandavg{1}.repetition);
sem_repetition_lowrt = std(grandavg{1}.repetition)/sqrt(length(subjects));

mean_repetition_highrt = mean(grandavg{2}.repetition);
sem_repetition_highrt = std(grandavg{2}.repetition)/sqrt(length(subjects));

clear grandavg

% Confidence 

for conf = 1:2

for sj = unique(subjects)
    
    % individual sj data
    data = alldata(alldata.subject == subjects(sj), :);
      
    nbins = 6;
      
    for u = 1:nbins
         
        if conf == 1 
            trls = find(data.prevDifficulty == u & data.prevConfidence <= 2);
        elseif conf == 2
            trls = find(data.prevDifficulty == u & data.prevConfidence >= 3);
        end
        
        % =========================================== %
        % history-dependent logistic
        % =========================================== %
        
        % previous response
        resps = [0 1];
        for r = 1:2
                        
            % take a further subset of trials with previous response 1 or 0
            thisdat         = data(trls,:);
            thisdat         = thisdat(thisdat.prevResponseDir == resps(r),:);
            
            [bias, slope, lapseLow, lapseHigh] = fitLogistic(...
                thisdat.currRelativeCoherence, thisdat.currResponseDir);
              
            % save betas
            grandavg{conf}.logisticHistory(sj, r, u, :) = [bias slope];
        end
               
    end
          
end

% ========================================================= %
% transform history-dependent bias into repetition bias
% ========================================================= %

resp1 = -1 * squeeze(grandavg{conf}.logisticHistory(:, 1, :, 1)); % take the 1st beta weigth only, ignore history-dependent slope
resp2 =      squeeze(grandavg{conf}.logisticHistory(:, 2, :, 1));

% since this is centred at 0.5, treat it that way
grandavg{conf}.repetition = (resp1 + resp2) ./ 2;

% ========================================================= %
% transform from log(odds) to probability
% ========================================================= %

logOdds2Probability = @(x) exp(x) ./ (1 + exp(x));
grandavg{conf}.repetition = logOdds2Probability(grandavg{conf}.repetition);

end

% save means and SEs for prev low and high conf

mean_repetition_lowconfidence = mean(grandavg{1}.repetition);
sem_repetition_lowconfidence = std(grandavg{1}.repetition)/sqrt(length(subjects));

mean_repetition_highconfidence = mean(grandavg{2}.repetition);
sem_repetition_highconfidence = std(grandavg{2}.repetition)/sqrt(length(subjects));

clear grandavg

% Correctness 

for corr = 1:2

for sj = unique(subjects)
    
    % individual sj data
    data = alldata(alldata.subject == subjects(sj), :);
      
    nbins = 6;
      
    for u = 1:nbins
         
        if corr == 1 
            trls = find(data.prevDifficulty == u & data.prevCorrect == 0);
        elseif corr == 2
            trls = find(data.prevDifficulty == u & data.prevCorrect == 1);
        end
        
        % =========================================== %
        % history-dependent logistic
        % =========================================== %
        
        % previous response
        resps = [0 1];
        for r = 1:2
                        
            % take a further subset of trials with previous response 1 or 0
            thisdat         = data(trls,:);
            thisdat         = thisdat(thisdat.prevResponseDir == resps(r),:);
            
            [bias, slope, lapseLow, lapseHigh] = fitLogistic(...
                thisdat.currRelativeCoherence, thisdat.currResponseDir);
              
            % save betas
            grandavg{corr}.logisticHistory(sj, r, u, :) = [bias slope];
        end
               
    end
          
end

% ========================================================= %
% transform history-dependent bias into repetition bias
% ========================================================= %

resp1 = -1 * squeeze(grandavg{corr}.logisticHistory(:, 1, :, 1)); % take the 1st beta weigth only, ignore history-dependent slope
resp2 =      squeeze(grandavg{corr}.logisticHistory(:, 2, :, 1));

% since this is centred at 0.5, treat it that way
grandavg{corr}.repetition = (resp1 + resp2) ./ 2;

% ========================================================= %
% transform from log(odds) to probability
% ========================================================= %

logOdds2Probability = @(x) exp(x) ./ (1 + exp(x));
grandavg{corr}.repetition = logOdds2Probability(grandavg{corr}.repetition);

end

% save means and SEs for prev incorrect and correct

mean_repetition_incorrect = mean(grandavg{1}.repetition);
sem_repetition_incorrect = std(grandavg{1}.repetition)/sqrt(length(subjects));

mean_repetition_correct = mean(grandavg{2}.repetition);
sem_repetition_correct = std(grandavg{2}.repetition)/sqrt(length(subjects));


% Correctness - P(repeat STIMULUS)

for corr = 1:2

for sj = unique(subjects)
    
    % individual sj data
    data = alldata(alldata.subject == subjects(sj), :);
      
    nbins = 6;
      
    for u = 1:nbins
         
        if corr == 1 
            trls = find(data.prevDifficulty == u & data.prevCorrect == 0);
        elseif corr == 2
            trls = find(data.prevDifficulty == u & data.prevCorrect == 1);
        end
        
        % =========================================== %
        % history-dependent logistic
        % =========================================== %
        
        % previous response
        resps = [-1 1];
        for r = 1:2
                        
            % take a further subset of trials with previous STIMULUS direction -1 or 1
            thisdat         = data(trls,:);
            thisdat         = thisdat(sign(thisdat.prevRelativeCoherence) == resps(r),:);
            
            [bias, slope, lapseLow, lapseHigh] = fitLogistic(...
                thisdat.currRelativeCoherence, thisdat.currResponseDir);
              
            % save betas
            grandavg{corr}.logisticHistory(sj, r, u, :) = [bias slope];
        end
               
    end
          
end

% ========================================================= %
% transform history-dependent bias into repetition bias
% ========================================================= %

resp1 = -1 * squeeze(grandavg{corr}.logisticHistory(:, 1, :, 1)); % take the 1st beta weigth only, ignore history-dependent slope
resp2 =      squeeze(grandavg{corr}.logisticHistory(:, 2, :, 1));

% since this is centred at 0.5, treat it that way
grandavg{corr}.repetition = (resp1 + resp2) ./ 2;

% ========================================================= %
% transform from log(odds) to probability
% ========================================================= %

logOdds2Probability = @(x) exp(x) ./ (1 + exp(x));
grandavg{corr}.repetition = logOdds2Probability(grandavg{corr}.repetition);

end

% save means and SEs for prev incorrect and correct

mean_repetition_stim_incorrect = mean(grandavg{1}.repetition);
sem_repetition_stim_incorrect = std(grandavg{1}.repetition)/sqrt(length(subjects));

mean_repetition_stim_correct = mean(grandavg{2}.repetition);
sem_repetition_stim_correct = std(grandavg{2}.repetition)/sqrt(length(subjects));

%% SIMULATED DATA %%

for iTeration = 1:iterations
    
    alldata.currResponseDir = simulatedResponses(:,iTeration);

    % RT

    for rt = 1:2

    for sj = unique(subjects)

        % individual sj data
        data = alldata(alldata.subject == subjects(sj), :);

        nbins = 6;

        for u = 1:nbins

            % first select trials for this bin
            trls = find(data.prevDifficulty == u);

            % find median rt
            datartsorted    = sortrows(data.prevRt(trls));
            medianrt        = datartsorted(floor(length(datartsorted)/2));

            % select low or high rts
            if rt == 1 
                trlssplit = trls(data.prevRt(trls) < medianrt);
            elseif rt == 2
                trlssplit = trls(data.prevRt(trls) >= medianrt);
            end

            % =========================================== %
            % history-dependent logistic
            % =========================================== %

            % previous response
            resps = [0 1];
            for r = 1:2

                % take a further subset of trials with previous response 1 or 0
                thisdat         = data(trlssplit,:);
                thisdat         = thisdat(thisdat.prevResponseDir == resps(r),:);

                [bias, slope, lapseLow, lapseHigh] = fitLogistic(...
                    thisdat.currRelativeCoherence, thisdat.currResponseDir);

                % save betas
                grandavg{rt}.logisticHistory(sj, r, u, :) = [bias slope];
            end    

        end

    end

    % ========================================================= %
    % transform history-dependent bias into repetition bias
    % ========================================================= %

    resp1 = -1 * squeeze(grandavg{rt}.logisticHistory(:, 1, :, 1)); % take the 1st beta weigth only, ignore history-dependent slope
    resp2 =      squeeze(grandavg{rt}.logisticHistory(:, 2, :, 1));

    % since this is centred at 0.5, treat it that way
    grandavg{rt}.repetition = (resp1 + resp2) ./ 2;

    % ========================================================= %
    % transform from log(odds) to probability
    % ========================================================= %

    logOdds2Probability = @(x) exp(x) ./ (1 + exp(x));
    grandavg{rt}.repetition = logOdds2Probability(grandavg{rt}.repetition);

    end

    % save means and SEs for prev low and high rt

    sim_mean_repetition_lowrt(iTeration,:) = mean(grandavg{1}.repetition);
    sim_sem_repetition_lowrt(iTeration,:) = std(grandavg{1}.repetition)/sqrt(length(subjects));

    sim_mean_repetition_highrt(iTeration,:) = mean(grandavg{2}.repetition);
    sim_sem_repetition_highrt(iTeration,:) = std(grandavg{2}.repetition)/sqrt(length(subjects));

    clear grandavg

    % Confidence 

    for conf = 1:2

     for sj = unique(subjects)

        % individual sj data
        data = alldata(alldata.subject == subjects(sj), :);

        nbins = 6;

        for u = 1:nbins
         
            if conf == 1 
                trls = find(data.prevDifficulty == u & data.prevConfidence <= 2);
            elseif conf == 2
                trls = find(data.prevDifficulty == u & data.prevConfidence >= 3);
            end

            % =========================================== %
            % history-dependent logistic
            % =========================================== %

            % previous response
            resps = [0 1];
            for r = 1:2

                % take a further subset of trials with previous response 1 or 0
                thisdat = data(trls,:);
                thisdat = thisdat(thisdat.prevResponseDir == resps(r),:);

                [bias, slope, lapseLow, lapseHigh] = fitLogistic(...
                    thisdat.currRelativeCoherence, thisdat.currResponseDir);

                % save betas
                grandavg{conf}.logisticHistory(sj, r, u, :) = [bias slope];
            end
               
        end
    
    end

    % ========================================================= %
    % transform history-dependent bias into repetition bias
    % ========================================================= %

    resp1 = -1 * squeeze(grandavg{conf}.logisticHistory(:, 1, :, 1)); % take the 1st beta weigth only, ignore history-dependent slope
    resp2 =      squeeze(grandavg{conf}.logisticHistory(:, 2, :, 1));

    % since this is centred at 0.5, treat it that way
    grandavg{conf}.repetition = (resp1 + resp2) ./ 2;

    % ========================================================= %
    % transform from log(odds) to probability
    % ========================================================= %

    logOdds2Probability = @(x) exp(x) ./ (1 + exp(x));
    grandavg{conf}.repetition = logOdds2Probability(grandavg{conf}.repetition);

    end

    % save means and SEs for prev low and high conf

    sim_mean_repetition_lowconfidence(iTeration,:) = mean(grandavg{1}.repetition);
    sim_sem_repetition_lowconfidence(iTeration,:) = std(grandavg{1}.repetition)/sqrt(length(subjects));

    sim_mean_repetition_highconfidence(iTeration,:) = mean(grandavg{2}.repetition);
    sim_sem_repetition_highconfidence(iTeration,:) = std(grandavg{2}.repetition)/sqrt(length(subjects));

    clear grandavg

    % Correctness 

    for corr = 1:2

      for sj = unique(subjects)

        % individual sj data
        data = alldata(alldata.subject == subjects(sj), :);

        nbins = 6;

        for u = 1:nbins

            if corr == 1 
                trls = find(data.prevDifficulty == u & data.prevCorrect == 0);
            elseif corr == 2
                trls = find(data.prevDifficulty == u & data.prevCorrect == 1);
            end

            % =========================================== %
            % history-dependent logistic
            % =========================================== %

            % previous response
            resps = [0 1];
            for r = 1:2

                % take a further subset of trials with previous response 1 or 0
                thisdat         = data(trls,:);
                thisdat         = thisdat(thisdat.prevResponseDir == resps(r),:);

                [bias, slope, lapseLow, lapseHigh] = fitLogistic(...
                    thisdat.currRelativeCoherence, thisdat.currResponseDir);

                % save betas
                grandavg{corr}.logisticHistory(sj, r, u, :) = [bias slope];
            end

        end

      end

    % ========================================================= %
    % transform history-dependent bias into repetition bias
    % ========================================================= %

    resp1 = -1 * squeeze(grandavg{corr}.logisticHistory(:, 1, :, 1)); % take the 1st beta weigth only, ignore history-dependent slope
    resp2 =      squeeze(grandavg{corr}.logisticHistory(:, 2, :, 1));

    % since this is centred at 0.5, treat it that way
    grandavg{corr}.repetition = (resp1 + resp2) ./ 2;

    % ========================================================= %
    % transform from log(odds) to probability
    % ========================================================= %

    logOdds2Probability = @(x) exp(x) ./ (1 + exp(x));
    grandavg{corr}.repetition = logOdds2Probability(grandavg{corr}.repetition);

    end

    % save means and SEs for prev low and high conf

    sim_mean_repetition_incorrect(iTeration,:) = mean(grandavg{1}.repetition);
    sim_sem_repetition_incorrect(iTeration,:) = std(grandavg{1}.repetition)/sqrt(length(subjects));

    sim_mean_repetition_correct(iTeration,:) = mean(grandavg{2}.repetition);
    sim_sem_repetition_correct(iTeration,:) = std(grandavg{2}.repetition)/sqrt(length(subjects));
    
    clear grandavg
    
    % Correctness - P(repeat STIMULUS)

    for corr = 1:2

      for sj = unique(subjects)

        % individual sj data
        data = alldata(alldata.subject == subjects(sj), :);

        nbins = 6;

        for u = 1:nbins

            if corr == 1 
                trls = find(data.prevDifficulty == u & data.prevCorrect == 0);
            elseif corr == 2
                trls = find(data.prevDifficulty == u & data.prevCorrect == 1);
            end

            % =========================================== %
            % history-dependent logistic
            % =========================================== %

            % previous response
            resps = [-1 1];
            for r = 1:2

                % take a further subset of trials with previous STIMULUS
                % DIRECTION -1 or 1
                thisdat         = data(trls,:);
                thisdat         = thisdat(sign(thisdat.prevRelativeCoherence) == resps(r),:);

                [bias, slope, lapseLow, lapseHigh] = fitLogistic(...
                    thisdat.currRelativeCoherence, thisdat.currResponseDir);

                % save betas
                grandavg{corr}.logisticHistory(sj, r, u, :) = [bias slope];
            end

        end

      end

    % ========================================================= %
    % transform history-dependent bias into repetition bias
    % ========================================================= %

    resp1 = -1 * squeeze(grandavg{corr}.logisticHistory(:, 1, :, 1)); % take the 1st beta weigth only, ignore history-dependent slope
    resp2 =      squeeze(grandavg{corr}.logisticHistory(:, 2, :, 1));

    % since this is centred at 0.5, treat it that way
    grandavg{corr}.repetition = (resp1 + resp2) ./ 2;

    % ========================================================= %
    % transform from log(odds) to probability
    % ========================================================= %

    logOdds2Probability = @(x) exp(x) ./ (1 + exp(x));
    grandavg{corr}.repetition = logOdds2Probability(grandavg{corr}.repetition);

    end

    % save means and SEs for prev low and high conf

    sim_mean_repetition_stim_incorrect(iTeration,:) = mean(grandavg{1}.repetition);
    sim_sem_repetition_stim_incorrect(iTeration,:) = std(grandavg{1}.repetition)/sqrt(length(subjects));

    sim_mean_repetition_stim_correct(iTeration,:) = mean(grandavg{2}.repetition);
    sim_sem_repetition_stim_correct(iTeration,:) = std(grandavg{2}.repetition)/sqrt(length(subjects));
    
    clear grandavg

end

%% Overall means and confidence intervals

% RT
sim_grand_mean_repetition_lowrt         = mean(sim_mean_repetition_lowrt);
    temp                                = sort(sim_mean_repetition_lowrt);
    sim_mean_repetition_lowrt_CI(1,:)   = temp(25,:);
    sim_mean_repetition_lowrt_CI(2,:)   = temp(975,:);
sim_grand_mean_repetition_highrt        = mean(sim_mean_repetition_highrt);
    temp                                = sort(sim_mean_repetition_highrt);
    sim_mean_repetition_highrt_CI(1,:)  = temp(25,:);
    sim_mean_repetition_highrt_CI(2,:)  = temp(975,:);

% Confidence
sim_grand_mean_repetition_lowconfidence         = mean(sim_mean_repetition_lowconfidence);
    temp                                        = sort(sim_mean_repetition_lowconfidence);
    sim_mean_repetition_lowconfidence_CI(1,:)   = temp(25,:);
    sim_mean_repetition_lowconfidence_CI(2,:)   = temp(975,:);
sim_grand_mean_repetition_highconfidence        = mean(sim_mean_repetition_highconfidence);
    temp                                        = sort(sim_mean_repetition_highconfidence);
    sim_mean_repetition_highconfidence_CI(1,:)  = temp(25,:);
    sim_mean_repetition_highconfidence_CI(2,:)  = temp(975,:);
    
% Correctness
sim_grand_mean_repetition_incorrect         = mean(sim_mean_repetition_incorrect);
    temp                                    = sort(sim_mean_repetition_incorrect);
    sim_mean_repetition_incorrect_CI(1,:)   = temp(25,:);
    sim_mean_repetition_incorrect_CI(2,:)   = temp(975,:);
sim_grand_mean_repetition_correct           = mean(sim_mean_repetition_correct);
    temp                                    = sort(sim_mean_repetition_correct);
    sim_mean_repetition_correct_CI(1,:)     = temp(25,:);
    sim_mean_repetition_correct_CI(2,:)     = temp(975,:);

% Correctness - P(repeat STIMULUS)
sim_grand_mean_repetition_stim_incorrect        = mean(sim_mean_repetition_stim_incorrect);
    temp                                        = sort(sim_mean_repetition_stim_incorrect);
    sim_mean_repetition_stim_incorrect_CI(1,:)  = temp(25,:);
    sim_mean_repetition_stim_incorrect_CI(2,:)  = temp(975,:);
sim_grand_mean_repetition_stim_correct          = mean(sim_mean_repetition_stim_correct);
    temp                                        = sort(sim_mean_repetition_stim_correct);
    sim_mean_repetition_stim_correct_CI(1,:)  	= temp(25,:);
    sim_mean_repetition_stim_correct_CI(2,:)    = temp(975,:);


    
%% Plots with simulation grand means, ALL simulation means, and real data

figure; hold on;

% p(repeat) per previous stimulus level, split for prev RT  
subplot(2,2,1); hold on;
for iTeration = 1:iterations
    plot(sim_mean_repetition_lowrt(iTeration,:),'Color',[0.7 0.7 1])
    plot(sim_mean_repetition_highrt(iTeration,:),'Color',[1 0.7 0.7])
end
v = plot(mean_repetition_lowrt,'b'); 
    errorbar(mean_repetition_lowrt,sem_repetition_lowrt,'.b');
w = plot(mean_repetition_highrt,'r');
    errorbar(mean_repetition_highrt,sem_repetition_highrt,'.r');
ylim([.3 .7]); 
xlim([0.25 6.75]); 
set(gca,'XTick',1:7,'XTickLabel',{'weakest','','','','','strongest'});
xlabel('Previous |evidence| strength')
ylabel('P(repeat)')
legend([v,w],'prev. low rt','prev. high rt','Location','southwest')
legend boxoff

% p(repeat) per previous stimulus level, split for prev confidence  
subplot(2,2,2); hold on;
for iTeration = 1:iterations
    plot(sim_mean_repetition_highconfidence(iTeration,:),'Color',[0.7 0.7 1])
    plot(sim_mean_repetition_lowconfidence(iTeration,:),'Color',[1 0.7 0.7])
end
v = plot(mean_repetition_highconfidence,'b'); 
    errorbar(mean_repetition_highconfidence,sem_repetition_highconfidence,'.b');
w = plot(mean_repetition_lowconfidence,'r');
    errorbar(mean_repetition_lowconfidence,sem_repetition_lowconfidence,'.r');
ylim([.3 .7]);
xlim([0.25 6.75]); 
set(gca,'XTick',1:7,'XTickLabel',{'weakest','','','','','strongest'});
ylabel('P(repeat)')
xlabel('Previous |evidence| strength')
legend([v,w],'prev. high confidence (3 or 4)', 'prev. low confidence (1 or 2)','Location','southwest')
legend boxoff;

% p(repeat) per previous stimulus level, split for prev correctness 
subplot(2,2,3); hold on;
for iTeration = 1:iterations
    plot(sim_mean_repetition_correct(iTeration,:),'Color',[0.7 0.7 1])
    plot(sim_mean_repetition_incorrect(iTeration,:),'Color',[1 0.7 0.7])
end
v = plot(mean_repetition_correct,'b'); 
    errorbar(mean_repetition_correct,sem_repetition_correct,'.b');
w = plot(mean_repetition_incorrect,'r');
    errorbar(mean_repetition_incorrect,sem_repetition_incorrect,'.r');
ylim([.3 .7]);
xlim([0.25 6.75]); 
set(gca,'XTick',1:7,'XTickLabel',{'weakest','','','','','strongest'});
ylabel('P(repeat)')
xlabel('Previous |evidence| strength')
legend([v,w],'prev. correct', 'prev. incorrect','Location','southwest')
legend boxoff;

% p(repeat STIMULUS) per previous stimulus level, split for prev correctness 
subplot(2,2,4); hold on;
for iTeration = 1:iterations
    plot(sim_mean_repetition_stim_correct(iTeration,:),'Color',[0.7 0.7 1])
    plot(sim_mean_repetition_stim_incorrect(iTeration,:),'Color',[1 0.7 0.7])
end
v = plot(mean_repetition_stim_correct,'b'); 
    errorbar(mean_repetition_stim_correct,sem_repetition_stim_correct,'.b');
w = plot(mean_repetition_stim_incorrect,'r');
    errorbar(mean_repetition_stim_incorrect,sem_repetition_stim_incorrect,'.r');
ylim([.3 .7]);
xlim([0.25 6.75]); 
set(gca,'XTick',1:7,'XTickLabel',{'weakest','','','','','strongest'});
ylabel('P(repeat STIMULUS)')
xlabel('Previous |evidence| strength')
legend([v,w],'prev. correct', 'prev. incorrect','Location','southwest')
legend boxoff;

set(gcf, 'Position',  [0, 0, 800, 800])

saveas(gcf,['P:/3018029.09/Graphs/p_repeat/' folder_str '/Compared_p(repeat) per previous stimulus level all lines.png'])
saveas(gcf,['P:/3018029.09/Graphs/p_repeat/' folder_str '/Compared_p(repeat) per previous stimulus level all lines.eps'])

    
%% Plots with simulation grand means, CIs, plus real data

figure; hold on;

% p(repeat) per previous stimulus level, split for prev RT  
subplot(2,2,1); hold on;
plot(sim_grand_mean_repetition_lowrt,'Color',[0.7 0.7 1])
    plot(sim_mean_repetition_lowrt_CI(1,:),'--','Color',[0.7 0.7 1])
    plot(sim_mean_repetition_lowrt_CI(2,:),'--','Color',[0.7 0.7 1])
plot(sim_grand_mean_repetition_highrt,'Color',[1 0.7 0.7])
    plot(sim_mean_repetition_highrt_CI(1,:),'--','Color',[1 0.7 0.7])
    plot(sim_mean_repetition_highrt_CI(2,:),'--','Color',[1 0.7 0.7])
x = plot(mean_repetition_lowrt,'b'); 
    errorbar(mean_repetition_lowrt,sem_repetition_lowrt,'.b');
y = plot(mean_repetition_highrt,'r');
    errorbar(mean_repetition_highrt,sem_repetition_highrt,'.r');
ylim([.3 .7]); 
xlim([0.25 6.75]); 
set(gca,'XTick',1:7,'XTickLabel',{'weakest','','','','','strongest'});
xlabel('Previous |evidence| strength')
ylabel('p(repeat)')
legend([x,y],'prev. low rt','prev. high rt','Location','southwest')
legend boxoff;

% p(repeat) per previous stimulus level, split for prev confidence  
subplot(2,2,2); hold on;
plot(sim_grand_mean_repetition_highconfidence,'Color',[0.7 0.7 1])
    plot(sim_mean_repetition_highconfidence_CI(1,:),'--','Color',[0.7 0.7 1])
    plot(sim_mean_repetition_highconfidence_CI(2,:),'--','Color',[0.7 0.7 1])
plot(sim_grand_mean_repetition_lowconfidence,'Color',[1 0.7 0.7])
    plot(sim_mean_repetition_lowconfidence_CI(1,:),'--','Color',[1 0.7 0.7])
    plot(sim_mean_repetition_lowconfidence_CI(2,:),'--','Color',[1 0.7 0.7])
x = plot(mean_repetition_highconfidence,'b'); 
    errorbar(mean_repetition_highconfidence,sem_repetition_highconfidence,'.b');
y = plot(mean_repetition_lowconfidence,'r'); 
    errorbar(mean_repetition_lowconfidence,sem_repetition_lowconfidence,'.r');
ylim([.3 .7]);
xlim([0.25 6.75]); 
set(gca,'XTick',1:7,'XTickLabel',{'weakest','','','','','strongest'});
ylabel('P(repeat)')
xlabel('Previous |evidence| strength')
legend([x,y],'prev. high confidence (3 or 4)', 'prev. low confidence (1 or 2)','Location','southwest')
legend boxoff;

% p(repeat) per previous stimulus level, split for prev correctness 
subplot(2,2,3); hold on;
plot(sim_grand_mean_repetition_correct,'Color',[0.7 0.7 1])
    plot(sim_mean_repetition_correct_CI(1,:),'--','Color',[0.7 0.7 1])
    plot(sim_mean_repetition_correct_CI(2,:),'--','Color',[0.7 0.7 1])
plot(sim_grand_mean_repetition_incorrect,'Color',[1 0.7 0.7])
    plot(sim_mean_repetition_incorrect_CI(1,:),'--','Color',[1 0.7 0.7])
    plot(sim_mean_repetition_incorrect_CI(2,:),'--','Color',[1 0.7 0.7])
x = plot(mean_repetition_correct,'b'); 
    errorbar(mean_repetition_correct,sem_repetition_correct,'.b');
y = plot(mean_repetition_incorrect,'r'); 
    errorbar(mean_repetition_incorrect,sem_repetition_incorrect,'.r');
ylim([.3 .7]);
xlim([0.25 6.75]); 
set(gca,'XTick',1:7,'XTickLabel',{'weakest','','','','','strongest'});
ylabel('P(repeat)')
xlabel('Previous |evidence| strength')
legend([x,y],'prev. correct', 'prev. incorrect','Location','southwest')
legend boxoff;

% p(repeat STIMULUS) per previous stimulus level, split for prev correctness 
subplot(2,2,4); hold on;
plot(sim_grand_mean_repetition_stim_correct,'Color',[0.7 0.7 1])
    plot(sim_mean_repetition_stim_correct_CI(1,:),'--','Color',[0.7 0.7 1])
    plot(sim_mean_repetition_stim_correct_CI(2,:),'--','Color',[0.7 0.7 1])
plot(sim_grand_mean_repetition_stim_incorrect,'Color',[1 0.7 0.7])
    plot(sim_mean_repetition_stim_incorrect_CI(1,:),'--','Color',[1 0.7 0.7])
    plot(sim_mean_repetition_stim_incorrect_CI(2,:),'--','Color',[1 0.7 0.7])
x = plot(mean_repetition_stim_correct,'b'); 
    errorbar(mean_repetition_stim_correct,sem_repetition_stim_correct,'.b');
y = plot(mean_repetition_stim_incorrect,'r'); 
    errorbar(mean_repetition_stim_incorrect,sem_repetition_stim_incorrect,'.r');
ylim([.3 .7]);
xlim([0.25 6.75]); 
set(gca,'XTick',1:7,'XTickLabel',{'weakest','','','','','strongest'});
ylabel('P(repeat STIMULUS)')
xlabel('Previous |evidence| strength')
legend([x,y],'prev. correct', 'prev. incorrect','Location','southwest')
legend boxoff;

set(gcf, 'Position',  [0, 0, 800, 800])

saveas(gcf,['P:/3018029.09/Graphs/p_repeat/' folder_str '/Compared_p(repeat) per previous stimulus level CIs.png'])
saveas(gcf,['P:/3018029.09/Graphs/p_repeat/' folder_str '/Compared_p(repeat) per previous stimulus level CIs.eps'])

  
     