1## Copyright (C) 2014 Tony Richardson
2##
3## This program is free software; you can redistribute it and/or modify it under
4## the terms of the GNU General Public License as published by the Free Software
5## Foundation; either version 3 of the License, or (at your option) any later
6## version.
7##
8## This program is distributed in the hope that it will be useful, but WITHOUT
9## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
11## details.
12##
13## You should have received a copy of the GNU General Public License along with
14## this program; if not, see <http://www.gnu.org/licenses/>.
15
16## -*- texinfo -*-
17## @deftypefn  {Function File} {[@var{h}, @var{pval}, @var{ci}, @var{stats}] =} ttest2 (@var{x}, @var{y})
18## @deftypefnx {Function File} {[@var{h}, @var{pval}, @var{ci}, @var{stats}] =} ttest2 (@var{x}, @var{y}, @var{Name}, @var{Value})
19## Test for mean of a normal sample with known variance.
20##
21## Perform a T-test of the null hypothesis @code{mean (@var{x}) ==
22## @var{m}} for a sample @var{x} from a normal distribution with unknown
23## mean and unknown std deviation.  Under the null, the test statistic
24## @var{t} has a Student's t distribution.
25##
26## If the second argument @var{y} is a vector, a paired-t test of the
27## hypothesis @code{mean (@var{x}) = mean (@var{y})} is performed.
28##
29## The argument @qcode{"alpha"} can be used to specify the significance level
30## of the test (the default value is 0.05).  The string
31## argument @qcode{"tail"}, can be used to select the desired alternative
32## hypotheses.  If @qcode{"alt"} is @qcode{"both"} (default) the null is
33## tested against the two-sided alternative @code{mean (@var{x}) != @var{m}}.
34## If @qcode{"alt"} is @qcode{"right"} the one-sided
35## alternative @code{mean (@var{x}) > @var{m}} is considered.
36## Similarly for @qcode{"left"}, the one-sided alternative @code{mean
37## (@var{x}) < @var{m}} is considered.  When @qcode{"vartype"} is @qcode{"equal"}
38## the variances are assumed to be equal (this is the default).  When
39## @qcode{"vartype"} is @qcode{"unequal"} the variances are not assumed equal.
40## When argument @var{x} is a matrix the @qcode{"dim"} argument can be
41## used to selection the dimension over which to perform the test.
42## (The default is the first non-singleton dimension.)
43##
44## If @var{h} is 0 the null hypothesis is accepted, if it is 1 the null
45## hypothesis is rejected. The p-value of the test is returned in @var{pval}.
46## A 100(1-alpha)% confidence interval is returned in @var{ci}. @var{stats}
47## is a structure containing the value of the test statistic (@var{tstat}),
48## the degrees of freedom (@var{df}) and the sample standard deviation
49## (@var{sd}).
50##
51## @end deftypefn
52
53## Author: Tony Richardson <richardson.tony@gmail.com>
54
55function [h, p, ci, stats] = ttest2(x, y, varargin)
56
57  alpha = 0.05;
58  tail  = 'both';
59  vartype = 'equal';
60
61  % Find the first non-singleton dimension of x
62  dim = min(find(size(x)~=1));
63  if isempty(dim), dim = 1; end
64
65  i = 1;
66  while ( i <= length(varargin) )
67    switch lower(varargin{i})
68      case 'alpha'
69        i = i + 1;
70        alpha = varargin{i};
71      case 'tail'
72        i = i + 1;
73        tail = varargin{i};
74      case 'vartype'
75        i = i + 1;
76        vartype = varargin{i};
77      case 'dim'
78        i = i + 1;
79        dim = varargin{i};
80      otherwise
81        error('Invalid Name argument.',[]);
82    end
83    i = i + 1;
84  end
85
86  if ~isa(tail, 'char')
87    error('Tail argument to ttest2 must be a string\n',[]);
88  end
89
90  m = size(x, dim);
91  n = size(y, dim);
92  x_bar = mean(x,dim)-mean(y,dim);
93  s1_var = var(x, 0, dim);
94  s2_var = var(y, 0, dim);
95
96  switch lower(vartype)
97    case 'equal'
98      stats.tstat = 0;
99      stats.df = (m + n - 2)*ones(size(x_bar));
100      sp_var = ((m-1)*s1_var + (n-1)*s2_var)./stats.df;
101      stats.sd = sqrt(sp_var);
102      x_bar_std = sqrt(sp_var*(1/m+1/n));
103    case 'unequal'
104      stats.tstat = 0;
105      se1 = sqrt(s1_var/m);
106      se2 = sqrt(s2_var/n);
107      sp_var = s1_var/m + s2_var/n;
108      stats.df = ((se1.^2+se2.^2).^2 ./ (se1.^4/(m-1) + se2.^4/(n-1)));
109      stats.sd = [sqrt(s1_var); sqrt(s2_var)];
110      x_bar_std = sqrt(sp_var);
111    otherwise
112      error('Invalid fifth (vartype) argument to ttest2\n',[]);
113  end
114
115  stats.tstat = x_bar./x_bar_std;
116
117  % Based on the "tail" argument determine the P-value, the critical values,
118  % and the confidence interval.
119  switch lower(tail)
120    case 'both'
121      p = 2*(1 - tcdf(abs(stats.tstat),stats.df));
122      tcrit = -tinv(alpha/2,stats.df);
123      %ci = [x_bar-tcrit*stats.sd; x_bar+tcrit*stats.sd];
124      ci = [x_bar-tcrit.*x_bar_std; x_bar+tcrit.*x_bar_std];
125    case 'left'
126      p = tcdf(stats.tstat,stats.df);
127      tcrit = -tinv(alpha,stats.df);
128      ci = [-inf*ones(size(x_bar)); x_bar+tcrit.*x_bar_std];
129    case 'right'
130      p = 1 - tcdf(stats.tstat,stats.df);
131      tcrit = -tinv(alpha,stats.df);
132      ci = [x_bar-tcrit.*x_bar_std; inf*ones(size(x_bar))];
133    otherwise
134      error('Invalid fourth (tail) argument to ttest2\n',[]);
135  end
136
137  % Reshape the ci array to match MATLAB shaping
138  if and(isscalar(x_bar), dim==2)
139    ci = ci(:)';
140    stats.sd = stats.sd(:)';
141  elseif size(x_bar,2)<size(x_bar,1)
142    ci = reshape(ci(:),length(x_bar),2);
143    stats.sd = reshape(stats.sd(:),length(x_bar),2);
144  end
145
146  % Determine the test outcome
147  % MATLAB returns this a double instead of a logical array
148  h = double(p < alpha);
149end
150