## Copyright (C) 2016 Andreas Stahel
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see .
## -*- texinfo -*-
## @deftypefn {Function File} {[@var{h}, @var{pval}, @var{ci}] =} binotest (@var{pos},@var{N},@var{p0})
## @deftypefnx {Function File} {[@var{h}, @var{pval}, @var{ci}] =} binotest (@var{pos},@var{N},@var{p0},@var{Name},@var{Value})
## Test for probability @var{p} of a binomial sample
##
## Perform a test of the null hypothesis @var{p} == @var{p0} for a sample
## of size @var{N} with @var{pos} positive results
##
##
## Name-Value pair arguments can be used to set various options.
## @qcode{"alpha"} can be used to specify the significance level
## of the test (the default value is 0.05). @qcode{"tail"}, can be used
## to select the desired alternative hypotheses. If the value is
## @qcode{"both"} (default) the null is tested against the two-sided
## alternative @code{@var{p} != @var{p0}}. If it is @qcode{"right"}
## the one-sided alternative @code{@var{p} > @var{p0}} is considered.
## Similarly for @qcode{"left"}, the one-sided alternative
## @code{@var{p} < @var{p0}} is considered.
##
## If @var{h} is 0 the null hypothesis is accepted, if it is 1 the null
## hypothesis is rejected. The p-value of the test is returned in @var{pval}.
## A 100(1-alpha)% confidence interval is returned in @var{ci}.
##
## @end deftypefn
## Author: Andreas Stahel
function [h, p, ci] = binotest(pos,n,p0,varargin)
% Set default arguments
alpha = 0.05;
tail = 'both';
i = 1;
while ( i <= length(varargin) )
switch lower(varargin{i})
case 'alpha'
i = i + 1;
alpha = varargin{i};
case 'tail'
i = i + 1;
tail = varargin{i};
otherwise
error('Invalid Name argument.',[]);
end
i = i + 1;
end
if ~isa(tail, 'char')
error('tail argument to vartest must be a string\n',[]);
end
% Based on the "tail" argument determine the P-value, the critical values,
% and the confidence interval.
switch lower(tail)
case 'both'
A_low = binoinv(alpha/2,n,p0)/n;
A_high = binoinv(1-alpha/2,n,p0)/n;
p = 2*min([binocdf(pos,n,p0), 1-binocdf(pos-1,n,p0)]);
p_low = fzero(@(pl)1-binocdf(pos-1,n,pl)-alpha/2,[0 1]);
p_high = fzero(@(ph) binocdf(pos,n,ph) -alpha/2,[0,1]);
ci = [p_low,p_high];
case 'left'
p = 1-binocdf(pos-1,n,p0);
p_high = fzero(@(ph) binocdf(pos,n,ph) -alpha,[0 1]);
ci = [0, p_high];
case 'right'
p = binocdf(pos,n,p0);
p_low = fzero(@(pl)1-binocdf(pos-1,n,pl) -alpha,[0 1]);
ci = [p_low 1];
otherwise
error('Invalid fifth (tail) argument to binotest\n',[]);
end
% Determine the test outcome
% MATLAB returns this a double instead of a logical array
h = double(p < alpha);
end
%!demo
%! % flip a coin 1000 times, showing 475 heads
%! % Hypothesis: coin is fair, i.e. p=1/2
%! [h,p_val,ci] = binotest(475,1000,0.5)
%! % Result: h = 0 : null hypothesis not rejected, coin could be fair
%! % P value 0.12, i.e. hypothesis not rejected for alpha up to 12%
%! % 0.0444 <= p <= 0.506 with 95% confidence
%!demo
%! % flip a coin 100 times, showing 65 heads
%! % Hypothesis: coin shows less than 50% heads, i.e. p<=1/2
%! [h,p_val,ci] = binotest(65,100,0.5,'tail','left','alpha',0.01)
%! % Result: h = 1 : null hypothesis is rejected, i.e. coin shows more heads than tails
%! % P value 0.0018, i.e. hypothesis not rejected for alpha up to 0.18%
%! % 0 <= p <= 0.76 with 99% confidence