% Response threshold model (Chapter 4, Box C)
% Written by David Sumpter

clear

n=40;
Paccept=0.1;   % Spontaneous decision to accept.
Quality1=1;     % Quality of sites
Quality2=0.5;   % Quality of sites

Pmax=1;    % Asymptotic probability of making a decision per time step.
              % 0<Pmax<1.
Pfind1=0.02;
Pfind2=0.02;


endtime=500;

T=10;      % Threshold response level.
k=1.0;     % Steepness of response

%Response function
figure(10)

a=[0:0.1:25];
f=Quality1*(Paccept+(Pmax-Paccept)*(a.^k./(T^k +a.^k)));
g=Quality2*(Paccept+(Pmax-Paccept)*(a.^k./(T^k +a.^k)));
plot(a,f,'r',a,g,'b')
axis([0 25 0 1])
xlabel('Number at site')
ylabel('Probability join')


%Number of trials stored in final distribution.
numtrials=100;

kvals=[1 9 0];


for i=1:numtrials

count=0;
for k=kvals
count=count+1;

e(1)=n;       %Number exploring
x1(1)=0;      %Stimulus left.
x2(1)=0;      %Stimulus right. 

t=1;
while ((e(t)>0) & (t<4000))
    
    %Set probabilities for 1 and 2
    P1 = Pfind1*Quality1*(Paccept+(Pmax-Paccept)*(x1(t)^k/(T^k +x1(t)^k)));
    P2 = Pfind2*Quality2*(Paccept+(Pmax-Paccept)*(x2(t)^k/(T^k +x2(t)^k)));
       
    if (k==0)
        P1 = Pfind1*Quality1;
        P2 = Pfind2*Quality2;
    end
    
    %Decide if they find and join a site for each exploring individual.
    e(t+1)=e(t);
    x1(t+1)=x1(t);
    x2(t+1)=x2(t);
    for j=1:e(t)
        r=rand;
        if (r<P1)
            x1(t+1)=x1(t+1)+1;
            e(t+1)=e(t+1)-1;
        elseif (r<P1+P2)
            x2(t+1)=x2(t+1)+1;
            e(t+1)=e(t+1)-1;
        end
    end
    t=t+1;    
end    

%Store end distributions
dx1(count,i)=x1(t);
dx2(count,i)=x2(t);
dt(count,i)=t;


if (i<=20)
    x1=[x1(1:t) ones(1000-t,1)'*x1(t)];
    x2=[x2(1:t) ones(1000-t,1)'*x2(t)];        
    sx1(count,i,:)=x1;
    sx2(count,i,:)=x2;
end

end

end

accuracy=shiftdim(mean(dx1./(dx1+dx2),4))
speed=shiftdim(mean(1./dt,4))


%Plot sequence of decisions and distribution of outcomes
figure(9)
    subplot(2,2,1)
    hold off
    plot(shiftdim(sx1(1,5,:),1),'k')
    hold on
    plot(shiftdim(sx2(1,5,:),1),'k:')
    ylabel('Number chosen option')
    xlabel('Time')
    title('(a)')
    %legend('Option 1','Option 2')
    axis([0 500 0 41])
subplot(2,2,2)
    hold off
    plot(shiftdim(sx1(2,2,:),1),'k')
    hold on
    plot(shiftdim(sx2(2,2,:),1),'k:')
    ylabel('Number chosen option')
    xlabel('Time')
    title('(b)')
    %legend('Option 1','Option 2')
    axis([0 500 0 41])
subplot(2,2,3)
[h x]=hist(dx1(1,:)/n,[0.05:0.05:0.95]);
bar(x,h/numtrials,'k')
xlabel('Proportion at X')
ylabel('frequency')
title('(c)')
axis([0 1 0 0.35])
subplot(2,2,4)
[h x]=hist(dx1(2,:)/n,[0.05:0.05:0.95]);
bar(x,h/numtrials,'k')
xlabel('Proportion at X')
ylabel('frequency')
title('(d)')
axis([0 1 0 0.35])

mean(dx1'/n)
mean(dt')
std(dt')

%Minimum time taken