%%
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
%%%%% Mockup DATA
%%% Stim lev. Probability of saying sound came from the left
dat=[
  -71.2870    2.0000   50.0000
  -61.1020    3.0000   50.0000
  -50.9190    3.0000   50.0000
  -40.7350    7.0000   50.0000
  -30.5510    7.0000   50.0000
  -20.3670   13.0000   50.0000
  -10.1840   25.0000   50.0000
         0   23.0000   50.0000
   10.1840   35.0000   50.0000
   20.3670   37.0000   50.0000
   30.5510   45.0000   50.0000
   40.7350   50.0000   50.0000
   50.9190   49.0000   50.0000
   61.1020   50.0000   50.0000
   71.2870   48.0000   50.0000];

datglm=dat;

%%  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)
plot(dat(:,1),(datglm(:,2)./datglm(:,3))*100,'x','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')
%%  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);

figure(2)
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 
%%
figure(4)
w=datglm(:,1);
poor=datglm(:,2);
total=datglm(:,3);

plot(w,poor./total,'x')
hold on
pad j = (poor+.5) ./ (total+1); 
plot(w,log(padj./(1-padj)),'x')




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)]