% define GrpSize outside of program
% define YesInGrp outside of program
% define NRep (number of simulations of the survey) outside of program
PrRealQ=0.5; % probability that dime lands Heads
PrYesDecoy=0.5; % probability that penny lands Heads
p=YesInGrp/GrpSize; % actual fraction of yes in group

%Each column in the matrices of this section represents one person in the group.
%We may assume that the first YesInGrp columns are the yes people.
TrueYes=[ones(NRep,YesInGrp) zeros(NRep,GrpSize-YesInGrp)];
RealQ=(rand(NRep,GrpSize)<PrRealQ); %each row tells who answers real question in one simulation
YesDecoy=(rand(NRep,GrpSize)<PrYesDecoy); %each row tells who would answer yes to decoy question

TrueYesAnswersYes=TrueYes.*(RealQ+(~RealQ).*YesDecoy);
TrueYesAnswersNo=TrueYes.*(~RealQ).*(~YesDecoy);
TrueNoAnswersYes=(~TrueYes).*(~RealQ).*(YesDecoy);
TrueNoAnswersNo=(~TrueYes).*(RealQ+(~RealQ).*(~YesDecoy));

TYY=TrueYesAnswersYes; %convenient abbreviation for name
TYN=TrueYesAnswersNo;
TNY=TrueNoAnswersYes;
TNN=TrueNoAnswersNo;

Summary=[sum(sum(TYY)) sum(sum(TNY)); sum(sum(TYN)) sum(sum(TNN))]/NRep;
%First row of this 2x2 table is the average number of yes answers in NRep
%runnings of the simulation, first for the true yes group, then for the true no group.
%Second row of this 2x2 table is the average number of no answers in NRep
%runnings of the simulation, first for the true yes group, then for the
%true no group.
NYes=sum(Summary(1,:));
%This is the average number of yes answers in NRep runnings of the
%simulation.
Estimate=sum(TYY'+TNY')/GrpSize; % naive estimate of p: fraction of yeses (in EACH simulation)
         % the text suggests that you could propose a better estimate!
         % MATLAB note: TYY' is the transpose of TYY, and 'sum' means sum
         % the columns.  The columns of TYY' are just the rows of TYY, and
         % the rows of TYY have a 1 for each 'yes' in the simulation.
         % Similarly for TNY.  Together these are all the 'yes' responses
         % in the simulation.
Bias=sum(Estimate-p)/NRep %find the average difference of the estimate from p over the simulations