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