%%
restoredefaultpath; matlabrc 
clear all
clc
close all

disp('Startup from auditory toolbox');
cd 'D:\COLLEGE_2007\STATS_PSYCHO\MATLABs\'
addpath( ...
    genpath(fullfile('D:','\MATL_DIR\psignifit')), ...
    'D:\MATL_DIR\BINOMIAL', ...
    'D:\MATL_DIR\CIBOOT');
cd  'D:\COLLEGE_2007\STATS_PSYCHO\MATLABs\'


%%
% See the documentation psych_options.m for a list of available options, 
% and refer to psych_glossary.m  for definitions and explanations where necessary. 
% The tutorial script psych_tool_demo.m begins to break down the operation of pfit, 
% demonstrating the use of some of the toolbox functions.
% Refer to psignifit/contents.m for a full list of available functions (
%%%%%%%%%%%type helpwin psignifit to browse the toolbox from within Matlab). 
%%%%%%%%%%      helpwin psignifit

clear all
close all
clc

%%% Stim lev. Probability Correct     number of measurements  yes no
dat=[
   -7.1287    0.0333   50.0000    2.0000   48.0000
   -6.1102    0.0556   50.0000    3.0000   47.0000
   -5.0919    0.0465   50.0000    3.0000   47.0000
   -4.0735    0.1304   50.0000    7.0000   43.0000
   -3.0551    0.1379   50.0000    7.0000   43.0000
   -2.0367    0.2500   50.0000   13.0000   37.0000
   -1.0184    0.4898   50.0000   25.0000   25.0000
         0    0.4595   50.0000   23.0000   27.0000
    1.0184    0.6818   50.0000   35.0000   15.0000
    2.0367    0.7391   50.0000   37.0000   13.0000
    3.0551    0.8889   50.0000   45.0000    5.0000
    4.0735    1.0000   50.0000   50.0000         0
    5.0919    0.9744   50.0000   49.0000    1.0000
    6.1102    1.0000   50.0000   50.0000         0
    7.1287    0.9565   50.0000   48.0000    2.0000];

dat(:,1)=dat(:,1)*10;

yes=ceil(dat(:,2).*dat(:,3));
no=dat(:,3)-yes;
dat=dat(:,1:3)
dat(:,2)=ceil(dat(:,2).*dat(:,3))
datglm=dat
dat(:,2)=dat(:,2)./dat(:,3)


%%  Fitting psychometric data  using glmfit
[b1,d1,s1] = glmfit(datglm(:,1),[datglm(:,2) datglm(:,3)],'binomial')
x = -100:.01:100; 
%y = glmval(b1,x,'probit',s1);
y = glmval(b1,x,'logit',s1);

figure(1)
hold off
plot(dat(:,1),(datglm(:,2)./datglm(:,3))*100,'x',x,y*100,'r-','LineWidth',2)
hold on
ylabel('Probability hearing sound from the left side','FontName','Arial','FontSize',13)
xlabel('Interaural time difference in microseconds','FontName','Arial','FontSize',13)
set(gca,'FontName','Arial','FontSize',13,'LineWidth',2);
grid 

p75=max(y(y>.749 & y<.751));
th75=x(find(y==p75))
plot(th75,p75*100,'o')
p25=max(y(y>.249 & y<.251));
th25=x(find(y==p25))
plot(th25,p25*100,'o')
%%







z_score=norminv(dat(:,2))

 

%%%%%y value is given by normpdf(norminv(0.980)) = .0478

figure
subplot(121)
plot(dat(:,1),dat(:,2)/max(dat(:,2)),'-o')

subplot(122)
plot(dat(:,2),z_score,'-o')



%%


%%% initial parameter setting
alpha=-0.2036 ;   
beta=2.633  ;  
gamma=0.03054 ;  
lambda=0.01412; 

%pfit(dat,'shape', 'cumulative gaussian', 'n_intervals', 1,'conf',[0.05 0.95],'cuts',[0.25 0.5 0.75], 'sens',15,'runs', 100)
pfit(dat,'shape', 'cumulative gaussian', 'n_intervals', 1,'conf',[0.05 0.95],'cuts',[0.25 0.5 0.75],'runs', 100)
[s rSeed STAT]=pfit(dat,'shape', 'cumulative gaussian', 'n_intervals', 1,'conf',[0.05 0.95],'cuts',[0.25 0.5 0.75], 'sens',15,'runs', 100)
%%%pfit(dat,'shape', 'cumulative gaussian', 'n_intervals', 1,'conf',[0.05 0.95],'cuts',[0.25 0.5 0.75], 'sens',10,'runs', 5000)

alpha=s.params.est(1);
beta=s.params.est(2);
gamma=s.params.est(3);
lambda=s.params.est(4);

%%BINOMIAL ERROR BARS
[lo up]=BinomialErrorBars(no,yes);

figure(1)
X=-8:0.01:8;
P = gamma + (1-gamma-lambda) * normcdf(X,alpha,beta);

figure(100)
hold off
plot(X,P*100)
            hold on
            errorbar(dat(:,1),dat(:,2)*100,abs(dat(:,2)-lo)*100,abs(up-dat(:,2))*100,'bo')
            
            
            
            
            
            
            
            %P = gamma + (1-gamma-lambda) * normcdf(X,alpha,1.2);
            %plot(X,P*100,'r-x')


%%% END PLOT 



% [db2,dd2,ds2] = glmfit(dat(:,1),[dat(:,2) dat(:,3)],'binomial');
% 
% %b = glmfit(x,[y n],'binomial','link','probit');
% yfit = glmval(db2, dat(:,1),'probit','size', dat(:,3));
% plot(dat(:,1), dat(:,2)./dat(:,3),'o',dat(:,1),yfit./dat(:,3),'-','LineWidth',2)

% b = glmfit(x,[y n],'binomial');
% yfit = glmval(b, x,'probit','size', n);
% figure(32)
% plot(x, y,'o',x,yfit*100,'-','LineWidth',2)

dat=datglm;

[b1,d1,s1] = glmfit(dat(:,1),[dat(:,2) dat(:,3)],'binomial')

x = -10:.01:10; 
y = glmval(b1,x,'logit',s1);
plot(dat(:,1),(dat(:,2)./dat(:,3))*100,'x',x,y*100,'r-')

set(gca,'FontName','Arial','FontSize',13,'LineWidth',2);





s.params.est(2)
s.params.lims(:,2)
% 
s.thresholds.est
s.thresholds.lims
% 
s.stats.deviance.D 


[s.stats.deviance.D 1-chi2cdf(s.stats.deviance.D ,15-4)]
[d1 1-chi2cdf(d1 ,15-2)]



