%psych_fit_plot.m

%to test: cd /Users/churchland/Dropbox/BehavData; load
%test_data_percent_correct; test_data = test_data_percent_correct;
%[xout yout B diagnotic_stuff] =psych_fit_plot(test_data(1,:),test_data(2,:),test_data(3,:))

function [xout, yout, B, diag] = psych_fit_plot(stimulus_strengths, number_of_correct, number_of_trials);

% stimulus_strengths = [0 0.2 0.4 0.6 0.8 1]‘;
% number_of_correct = [0 2 25 48 50 48]‘;
% number_of_correct = [15 20 25 48 50 48]‘;
% number_of_trials = 50 .* ones(1,length(stimulus_intensities))’;

data = [stimulus_strengths' number_of_correct' number_of_trials'] ;

priors.m_or_a = ’None’;
priors.w_or_b = ’None’;
priors.lambda = ’Uniform(0,.1)’;
priors.gamma = ’Uniform(0,.1)’;

%B = BayesInference(data, priors, ’sigmoid’, ‘gauss’, ‘core’, ‘ab’, ‘nafc’, 1);
B = BootstrapInference(data, priors, ’sigmoid’, ‘gauss’, ‘core’, ‘ab’, ‘nafc’, 1);
%B = BootstrapInference(data,priors);

%This next little bit is just because the psignifit people changed the name of the fitted parameters from thetahat to params_estimate for the newer version;
if isfield(B,’thetahat’)
diag = Diagnostics(B.data, B.thetahat, ’sigmoid’, B.sigmoid, ’core’, B.core, ’nafc’, B.nafc, ’cuts’, B.cuts);
elseif isfield(B,’params_estimate’)
diag = Diagnostics(B.data, B.params_estimate, ’sigmoid’, B.sigmoid, ’core’, B.core, ’nafc’, B.nafc, ’cuts’, B.cuts); %for params_estimate(1) = mean, params_estimate(2) = SD, params_estimate(3:4) = deviation from 0 for left and right side of fits.
end;

% confidence intervals
% ci = getCI(B, 2, 0.95);
% the fit
% plot(diag.pmf(:,1), diag.pmf(:,2), ‘-’);
% hold on;
% the data points
% plot(B.data(:,1), B.data(:,2)./B.data(:,3), ‘o’);

xout = diag.pmf(:,1);
yout = diag.pmf(:,2);

plot_me= 1;
if plot_me == 1
%  subplot(1,2,1)
plot(stimulus_strengths, number_of_correct./number_of_trials,’ko’); hold on;
plot(diag.pmf(:,1),diag.pmf(:,2),’k-’)
axprefs(gca);ylabel(‘Percent correct’);
xlabel(‘stimulus strength’);
end;

end