## 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