octave_packages/m/statistics/tests/t_test_regression.m

   1 ## Copyright (C) 1995-2012 Kurt Hornik
   2 ##
   3 ## This file is part of Octave.
   4 ##
   5 ## Octave is free software; you can redistribute it and/or modify it
   6 ## under the terms of the GNU General Public License as published by
   7 ## the Free Software Foundation; either version 3 of the License, or (at
   8 ## your option) any later version.
   9 ##
  10 ## Octave is distributed in the hope that it will be useful, but
  11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
  12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 ## General Public License for more details.
  14 ##
  15 ## You should have received a copy of the GNU General Public License
  16 ## along with Octave; see the file COPYING.  If not, see
  17 ## <http://www.gnu.org/licenses/>.
  18
  19 ## -*- texinfo -*-
  20 ## @deftypefn {Function File} {[@var{pval}, @var{t}, @var{df}] =} t_test_regression (@var{y}, @var{x}, @var{rr}, @var{r}, @var{alt})
  21 ## Perform an t test for the null hypothesis @code{@var{rr} * @var{b} =
  22 ## @var{r}} in a classical normal regression model @code{@var{y} =
  23 ## @var{x} * @var{b} + @var{e}}.  Under the null, the test statistic @var{t}
  24 ## follows a @var{t} distribution with @var{df} degrees of freedom.
  25 ##
  26 ## If @var{r} is omitted, a value of 0 is assumed.
  27 ##
  28 ## With the optional argument string @var{alt}, the alternative of
  29 ## interest can be selected.  If @var{alt} is @code{"!="} or
  30 ## @code{"<>"}, the null is tested against the two-sided alternative
  31 ## @code{@var{rr} * @var{b} != @var{r}}.  If @var{alt} is @code{">"}, the
  32 ## one-sided alternative @code{@var{rr} * @var{b} > @var{r}} is used.
  33 ## Similarly for @var{"<"}, the one-sided alternative @code{@var{rr} *
  34 ## @var{b} < @var{r}} is used.  The default is the two-sided case.
  35 ##
  36 ## The p-value of the test is returned in @var{pval}.
  37 ##
  38 ## If no output argument is given, the p-value of the test is displayed.
  39 ## @end deftypefn
  40
  41 ## Author: KH <Kurt.Hornik@wu-wien.ac.at>
  42 ## Description: Test one linear hypothesis in linear regression model
  43
  44 function [pval, t, df] = t_test_regression (y, x, rr, r, alt)
  45
  46   if (nargin == 3)
  47     r   = 0;
  48     alt = "!=";
  49   elseif (nargin == 4)
  50     if (ischar (r))
  51       alt = r;
  52       r   = 0;
  53     else
  54       alt = "!=";
  55     endif
  56   elseif (! (nargin == 5))
  57     print_usage ();
  58   endif
  59
  60   if (! isscalar (r))
  61     error ("t_test_regression: R must be a scalar");
  62   elseif (! ischar (alt))
  63     error ("t_test_regression: ALT must be a string");
  64   endif
  65
  66   [T, k] = size (x);
  67   if (! (isvector (y) && (length (y) == T)))
  68     error ("t_test_regression: Y must be a vector of length rows (X)");
  69   endif
  70   s      = size (rr);
  71   if (! ((max (s) == k) && (min (s) == 1)))
  72     error ("t_test_regression: RR must be a vector of length columns (X)");
  73   endif
  74
  75   rr     = reshape (rr, 1, k);
  76   y      = reshape (y, T, 1);
  77   [b, v] = ols (y, x);
  78   df     = T - k;
  79   t      = (rr * b - r) / sqrt (v * rr * inv (x' * x) * rr');
  80   cdf    = tcdf (t, df);
  81
  82   if (strcmp (alt, "!=") || strcmp (alt, "<>"))
  83     pval = 2 * min (cdf, 1 - cdf);
  84   elseif strcmp (alt, ">")
  85     pval = 1 - cdf;
  86   elseif strcmp (alt, "<")
  87     pval = cdf;
  88   else
  89     error ("t_test_regression: the value `%s' for alt is not possible", alt);
  90   endif
  91
  92   if (nargout == 0)
  93     printf ("pval: %g\n", pval);
  94   endif
  95
  96 endfunction