function [ME,D,Dprob,halfCI]=pollsig(N,V,Trials,details)
% How significant are differences in poll results?
% format [ME,D,Dprob,halfCI]=pollsig(N,V,Trials,details);
% Input:
%  Required:
%    N = Number of individuals polled
%    V=  Column or row vector with proportion of responses, 
%        need not sum to 1, but sum(V)<1;
%        In this implementation only the two most common items will be 
%        tested.
%  Optional:
%    Trials=Number of Monte Carlo trials used to judge significance
%       if Trials not specified, 1e4 trials will be used
%       if Trials=0, then the covariance formula will be used to judge
%                    significance
%    details=0, supress all output within the m.file.
% Output:
%   ME=Margin of Error for the poll
%   D=difference in proportions, length(V) x length(V) symmetric matrix.
%   Dprob=two-sided p value for test of equality of proportions of the 2 %
%         most common items
%         Dprob will have a minimum value of 1/Trials
%   halfCI=half the 95% CI for difference in proportions of the two most
%   common items.
%  Reference: Larsen & Marx (2006) Introduction to Mathematical Statistics
%  4th Edition p. 372

% Written in 2003 for EEOS601 by Eugene.Gallagher@umb.edu
% Dept of Environmental, Earth & Ocean Sciences.

if nargin<4
    details=1;
    MC=1;
    if nargin<3
        Trials=1e4;
    elseif Trials==0
        MC=0;
    end
end

% Calculate Margin of Error, p. 372, Larsen & Marx (2006) 4th edition
ME=norminv(0.975)/(2*sqrt(N));

if details
    fprintf('The margin of error for a poll of size %d is %3.1f%%.\n',...
        N,ME*100);
end

% Monte Carlo simulation

if details;
    fprintf('\nMonte Carlo simulation based on %d trials:\n',Trials);
end
V=V(:); % Change V to a column vector
V=flipud(sort(V));V=V(1:2);  % This m.file will only calculate significance
                             % for top two categories.
tallys=zeros(Trials,2); % A column vector with rows=Trials & 2 columns; 
                        % tallys differences
tallyHo=zeros(Trials,2);% This will store the results for testing Ho:
                        % p1=p2;
ExpP=mean(V);
for i=1:Trials
    poll=rand(N,1);     % Creates a vector with uniformly distributed
                        % random numbers on the interval 0,1
    tallys(i,1)=sum(poll<=V(1));
    tallys(i,2)=sum( (poll>V(1)) & (poll <= (V(1)+ V(2))) );
    tallyHo(i,1)=sum(poll<=ExpP);
    tallyHo(i,2)=sum( (poll>ExpP) & (poll <= 2*ExpP));
end
DifferenceHo = (tallyHo(:,1) - tallyHo(:,2))/N; % Calculate the differences
% for all Trials under Ho: p1=p2
D=abs(V(1)-V(2));
Dprob=max([1 sum(abs(DifferenceHo)>=D)])/Trials;
if details & Dprob<0.001  % change the format so that it is in full form 
                          % only for low p values:
    fprintf(...
  'Under the null hypothesis of equal proportions and %d trials,\n',Trials)
    fprintf(...
 'the 2-sided prob. of observing a %5.3f%% difference by chance = %d\n',...
        D*100,Dprob);
elseif details
    fprintf(...
'Under the null hypothesis of equal proportions, and %d trials,\n',...
Trials)
    fprintf(...
'the 2-sided P of observing a %4.2f%% difference by chance = %5.3f\n',...
        D*100,Dprob);
end
Diff = (tallys(:,1) - tallys(:,2))/N;
% 95% CI via Monte Carlo simulation
sortedDiff=sort(Diff);
lMC95CIpi=floor(0.025*Trials); % find the index for the lower 95% cutoff
uMC95CIpi=ceil(0.975*Trials);  % find the index for the upper 95% cutoff
medpi=round(0.5*Trials);       % find the median, should be close to or 
                               % identical to the expected value.
% Save the three outputs in the row vector DLowExpUp
DLowExpUp=[sortedDiff(lMC95CIpi) sortedDiff(medpi) sortedDiff(uMC95CIpi)];
halfCI=(DLowExpUp(3)-DLowExpUp(1))/2;
if details
    fprintf(...
'\nLower 95%% confidence limit, median, and upper 95%% confidence limit based on %d trials:\n',Trials)
    fprintf(...
'Lower 95%% CI  \tMedian \tUpper 95%% CI \n \t%4.2f%% \t\t%4.2f%% \t%4.2f%% \n',...
       DLowExpUp(1)*100,DLowExpUp(2)*100,DLowExpUp(3)*100)
   fprintf('\nDifference +/- half 95%% CI:  %4.1f%% +/- %4.1f%%\n',D*100,halfCI*100) 
end