function [a, P]=winstay(EstQ,P,beta)
% function [a,P]=winstay(EstQ,P,beta)
% Input: EstQ, Estimated values of the payoffs
% P = Probabilities of playing each machine
% beta= parameter to adjust the probabilities, between 0 and 1
% Output: a = Which machine to play
% P = Probabilities for each machine
[vals,idx]=max(EstQ);
winner=idx(1); %Index of our "winning" machine
%Update the probabilities. We need to do P(winner) separately.
NumMachines=length(P);
P(winner)=P(winner)+beta*(1-P(winner));
Temp=1:NumMachines;
Temp(winner)=[]; %Temp now holds the indices of all "losers"
P(Temp)=(1-beta)*P(Temp);
%Probabilities are all updated- Choose machine a w/prob P(a)
x=rand;
TotalBins=histc(x,[0,cumsum(P)']);
a=find(TotalBins==1);