%Lecture_041112.m

clear all
load resultsStaircase

pInit.t = .1;
pInit.b = 2;
pInit.shutup = 'please';
freeList ={'t','b'};

[pBest,logLikelihoodBest] = fit('fitPsychometricFunction',pInit,freeList,results,'Weibull');

figure(1)
clf
plotPsycho(results,'Coherence',pBest,'Weibull')

%% generate a 'fake' data set based on the best fitting Weibul.

intensities = unique(results.intensity);
p = Weibull(pBest,intensities);

fakeResults.intensity = results.intensity;
fakeResults.response = floor(rand(size(results.response))+ Weibull(pBest,results.intensity));
[pFakeBest,logLikelihoodBest] = fit('fitPsychometricFunction',pInit,freeList,fakeResults,'Weibull');
figure(2)
clf
plotPsycho(fakeResults,'Coherence',pFakeBest,'Weibull')


figure(2)
%% do this 1000 times to get an estimate of the distribution of our threshold

nReps = 300;
threshList = zeros(1,nReps);

%set up the wait bar
h = waitbar(0,'Please wait...');

p = Weibull(pBest,results.intensity);

for i=1:nReps
    fakeResults.response = floor(rand(size(results.response))+p);
    [pFakeBest,logLikelihoodBest] = fit('fitPsychometricFunction',...
        pInit,freeList,fakeResults,'Weibull');
    threshList(i) = pFakeBest.t;
    
    waitbar(i/nReps,h)
    
end
close(h);

%% histogram of thresholds

figure(2)
clf
hist(threshList,50);
hold on
ylim = get(gca,'YLim');
plot(pBest.t*[1,1],ylim,'r-','LineWidth',2);


%% confidence interval
CIpct = 68;

CI = [prctile(threshList,(100-CIpct)/2),...
prctile(threshList,100-(100-CIpct)/2)];

figure(3)
clf
hold on
plot(1,log(pBest.t),'ko','MarkerFaceColor','b');
errorbar(1,log(pBest.t),log(CI(1)),log(CI(2)),'b');
set(gca,'YTick',log([.01,.05,.1,.5]));
logy2raw
set(gca,'XLim',[0,2]);












