## 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 <http://www.gnu.org/licenses/>.

## -*- 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 <Andreas.Stahel@bfh.ch.
## based on the code ttest.m by Tony Richardson <richardson.tony@gmail.com>

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

