% PA_HOW2RECIPROBIT

% 2013 Marc van Wanrooij
% e: marcvanwanrooij@neural-code.com

close all
clearvars
clc

% Check whether my custom functions are available
hasCbrewer  = ~isempty(which('cbrewer'));
hasNicegraph = ~isempty(which('nicegraph'));
hasSavegraph = ~isempty(which('savegraph'));


%% Histogram
fname	= fullfile('..','..','data','reactiontime.mat');
load(fname); % load the reaction time data into Matlab workspace
rt		= RT.easy; % let's take the reaction times (ms) for the "easy" task
x		= 0:50:1000; % reaction time bins (ms)
dx		= x(2) - x(1);
edges	= [x - dx/2, x(end) + dx/2];

N	= histcounts(rt,edges);  % count reactions in bin
N	= 100*N./sum(N); % normalization / percentage

figure(1) % graphics
h	= bar(x,N); % histogram/bar plot
% Making figure nicer to look at
set(h,'FaceColor',[.7 .7 .7]); % setting the color of the bars to gray
box off; % remove the "box" around the graph
axis square;
xlabel('reaction time (ms)'); % and setting labels is essential for anyone to understand graphs!
ylabel('probability');
xlim([0 1000]);
set(gca,'YTick',[]);
nicegraph;

mean(rt)
% optional, save figure
fname = fullfile('..','..','images',[mfilename '_histogram_ori.svg']);
savegraph(fname,'svg');

%%
[h,p]	= kstest(zscore(rt)) % for reaction time



%% Inverse reaction time
rtinv	= 1000./rt; % inverse reaction time / promptness (ms-1)

n		= numel(x); % number of bins in reaction time plot
x		= linspace(0,10,n); % promptness bins
dx		= x(2) - x(1);
edges	= [x - dx/2, x(end) + dx/2];

N	= histcounts(rtinv,edges);  % count reactions in bin

N		= N./sum(N);
figure(2)
clf
h = bar(x,N);
hold on

set(h,'FaceColor','#eeb4d7','EdgeColor','#d74ea2');
% box off
% axis square;
xlabel('promptness (s^{-1})');
ylabel('probability');
title('reciprocal time axis');
set(gca,'YTick',0:5:100,'XTick',0:1:8);
xlim([0 8]);

% Does this look like a Gaussian?
% Let's plot a Gaussian curve with mean and standard deviation from the
% promptness data
mu	= mean(rtinv);
sd	= std(rtinv);
y	= normpdf(x,mu,sd);
y	= y./sum(y);
xi		= linspace(0,10,1000); % promptness bins
yi	= normpdf(xi,mu,sd);
yi	= yi./sum(yi)*1000/n*1.05;
plot(xi,yi,'k-','LineWidth',2);

plot(x,y,'ks','LineWidth',2,'MarkerFaceColor','w','MarkerSize',15);

xlabel('promptness (s^{-1})'); % and setting labels is essential for anyone to understand graphs!
ylabel('probability');
xlim([0 9]);
set(gca,'YTick',[]);
nicegraph;

applyTheme('barbie');
% optional, save figure
fname = fullfile('..','..','images',[mfilename '_promptness_histogram_ori.svg']);
savegraph(fname,'svg');


% And test it with a one-sample kolmogorov-smirnov test
% Since this test compares with the normal distribution, we have to
% normalize the data, which we can do with zscore, which is the same as:
% z = (x-mean(x))/std(x)
[h,p]	= kstest(zscore(rt)); % for reaction time
if h
	str		= ['The null hypothesis that the reaction time distribution is Gaussian distributed is rejected, with P = ' num2str(p)];
else
	str		= ['The null hypothesis that the reaction time distribution is Gaussian distributed is NOT rejected, with P = ' num2str(p)];
end
disp(str); % display the results in command window
[h,p]	= kstest(zscore(rtinv)); % for inverse reaction time
if h
	str		= ['The null hypothesis that the inverse reaction time distribution is Gaussian distributed is rejected, with P = ' num2str(p)];
else
	str		= ['The null hypothesis that the inverse reaction time distribution is Gaussian distributed is NOT rejected, with P = ' num2str(p)];
end
disp(str); % display the results in command window



%% Cumulative probability
% Normalized scale and nicer shape
x = sort(rtinv);
n = numel(rtinv); % number of data points
y = (1:n)./n; % cumulative probability for every data point
figure(3)
subplot(121)
plot(x,y,'k.');
hold on

% Now, plot it cumulative probability in quantiles
% this is easier to compare between different distributions
p		= [1 2 5 10:10:90 95 98 99]/100; % some arbitrary probabilities
q		= quantile(rtinv,p); % instead of hist, let's use quantiles

h = plot(q,p,'ko','LineWidth',1,'MarkerFaceColor','r','MarkerSize',12);
hold on
xlabel('promptness (s^{-1})');
ylabel('cumulative probability');
title('cumulative probability of reciprocal reaction time plot');
box off
axis square;
set(gca,'YTick',0:0.1:1,'XTick',0:1:10);
yline(0.5);
xlim([0 10]);
legend(h,'Quantiles','Location','SE');
nicegraph;
bf_text(0.05,0.95,'a','FontSize',21);






% Probit
figure(3)
subplot(122)
cla
% raw data
x = -1000./sort((rt)); % multiply by -1 to mirror abscissa
n = numel(rtinv); % number of data points
y = probit((1:n)./n); % cumulative probability for every data point converted to probit scale
plot(x,y,'k.');
hold on

% quantiles
p		= [1 2 5 10:10:90 95 98 99]/100;
prbt	= probit(p);
q		= quantile(rt,p);
q		= -1000./q;
xtick	= sort(-1000./(oct2bw(100,0:7))); % some arbitrary xticks

plot(q,prbt,'ko','Color','k','MarkerFaceColor','r','LineWidth',2,'MarkerSize',12);
hold on

set(gca,'XTick',xtick,'XTickLabel',-1000./xtick);
xlim([min(xtick) max(xtick)]);
set(gca,'YTick',prbt,'YTickLabel',p);
ylim([probit(0.1/100) probit(99.9/100)]);

xlabel('reaction time (ms)');
ylabel('cumulative probability');
title('probit ordinate with reciprocal abscissa');
nicegraph;
% this should be a straight line
x = q;
y = prbt;
b = regstats(y,x);
h = regline(b.beta,'k-');
set(h,'Color','r','LineWidth',2);
bf_text(0.05,0.95,'b','FontSize',21)

% optional, save figure
fname = fullfile('..','..','images',[mfilename '_promptness_cdf_ori.svg']);
savegraph(fname,'svg');