%Lecture_040912.m

%% Demonstrate the Central Limit Theorem:

% Single sample:
n =25;  %sample size

x = randn(n,1);  %sample from standard normal (mean 0 s.d. 1)

figure(1)
clf
hist(x,50);
title(sprintf('mean = %5.2f',mean(x)));

%% draw a histogram of the sampling distribution of the mean

nReps =1000;

y = randn(n,nReps);
sampDist = mean(y);

figure(2)
clf
hist(sampDist,50);

title(sprintf('mean = %5.2f, s.d. = %5.2f',mean(sampDist),std(sampDist)));

% should have a standard deviation that is 1/sqrt(n);
%%  Function handle for our statistic

myStat = @mean;

%our one sample:

x = randn(n,1)+.25;


%% Sampling with replacement

% n x nReps matrix of values sampled from x
mySamp =   x(ceil(rand(n,nReps)*25));

% take the mean of the 1000 columns
ySamp = myStat(mySamp);

% show the means as a histogram.  compare to figure 1.  Same SD, but the
% mean of the distribution should be near the mean of the sample.

figure(3)
clf
hist(ySamp,50);

pBoot = sum(ySamp<0)/nReps;

title(sprintf('mean = %5.2f, s.d. = %5.2f, %5.2f%% are less than zero',mean(ySamp),std(ySamp),100*pBoot));



CI = prctile(ySamp,5);
figure(3)
hold on
ylim = get(gca,'YLim');
plot(CI*[1,1],ylim,'r-','LineWidth',2);

%% compare to t-test (if you have the stats toolbox)

 [H,P,CI,STATS] = ttest(x,0,.05,'right');

 title(sprintf('bootstrap p-value %5.4f, t-test p-value %5.4f',pBoot,P));


%% my 'bootstrap' function
figure(3)
[CI,sampleStat] = bootstrap(myStat,x,1000,90);
hold on
plot(CI(1)*[1,1],ylim,'g-','LineWidth',2);

%% Next lecture: bootrapping to estimate standard error for psychophysical thresholds.

load resultsStaircase

%find the acutal threshold

pInit.t = .2;
pInit.b = 2;

[pBest,logLikelihoodBest] = fit('fitPsychometricFunction',pInit,{'t','b'},results,'Weibull')
plotPsycho(results,'Coherence',pBest,'Weibull')

