%Homework_2_Boynton

%Choose 'internal' Weibull parameters:

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

%Choose a list of numbers of trials (per intensity)
nTrialList = [1,2,4,8,16,32,64,128];

%Choose number of repetitions for the simulation
nReps = 2000;

%Choose intensity values:
intensityList = [.02,.04,.08,.16,.32,.64];

prob = Weibull(p,intensityList);

%results.intensity = intensityList;

nTrials = nTrialList*length(intensityList);

CI = zeros(length(nTrialList),2);

for i = 1:length(nTrialList)
    disp(sprintf('Length: %3d (%d of %d)',nTrials(i),i,length(nTrialList)));
    tmp = repmat(intensityList,nTrialList(i),1);
    results.intensity= tmp(:)';
    prob = Weibull(p,results.intensity);
    results.response = prob;
    CI(i,:) = bootstrapWeibullThreshold(results,p,nReps);
end
%%
save data %to avoid running this 1/2 hour program again

%%
load data

% Calculate the half-width for each CI
halfWidth = ((CI(:,2)-CI(:,1))/2)';

% fit the (log of the) results with a linear function (first order polynomial)
poly = polyfit(log(nTrials),log(halfWidth),1);

% Evaluate the polynomial for plotting
y = polyval(poly,log(nTrials));

% Plot the data and polynomial fit
figure(1)
clf
plot(log(nTrials),log(halfWidth),'o','MarkerFaceColor','k')
hold on
plot(log(nTrials),y,'k-')
set(gca,'XTick',log(nTrials));
set(gca,'YTick',log([.001,.01,.02,.05,.1,.2]));
logx2raw
logy2raw
xlabel('Number of trials');
ylabel('Half-width of CI');

title(sprintf('Slope = %5.3f',poly(1)));