view scripts/statistics/tests/manova.m @ 19074:0850b5212619

doc: Add @nospell macro around proper names in documentation. * aspell-octave.en.pws: Remove proper names from Octave's custom dictionary. * doc/interpreter/audio.txi, doc/interpreter/diagperm.txi, doc/interpreter/diffeq.txi, doc/interpreter/emacs.txi, doc/interpreter/func.txi, doc/interpreter/geometry.txi, doc/interpreter/intro.txi, doc/interpreter/linalg.txi, doc/interpreter/numbers.txi, doc/interpreter/preface.txi, doc/interpreter/quad.txi, doc/interpreter/sparse.txi, libinterp/corefcn/colloc.cc, libinterp/corefcn/ellipj.cc, libinterp/corefcn/gammainc.cc, libinterp/corefcn/hess.cc, libinterp/corefcn/quadcc.cc, libinterp/corefcn/rand.cc, libinterp/corefcn/sqrtm.cc, libinterp/dldfcn/ccolamd.cc, libinterp/dldfcn/colamd.cc, libinterp/dldfcn/dmperm.cc, libinterp/dldfcn/symrcm.cc, scripts/general/quadgk.m, scripts/general/quadl.m, scripts/linear-algebra/commutation_matrix.m, scripts/linear-algebra/condest.m, scripts/linear-algebra/duplication_matrix.m, scripts/linear-algebra/expm.m, scripts/linear-algebra/housh.m, scripts/linear-algebra/krylov.m, scripts/linear-algebra/onenormest.m, scripts/linear-algebra/qzhess.m, scripts/optimization/fminsearch.m, scripts/optimization/glpk.m, scripts/signal/arch_fit.m, scripts/signal/bartlett.m, scripts/signal/blackman.m, scripts/signal/diffpara.m, scripts/signal/durbinlevinson.m, scripts/signal/hamming.m, scripts/signal/hanning.m, scripts/sparse/eigs.m, scripts/sparse/pcg.m, scripts/sparse/pcr.m, scripts/specfun/ellipke.m, scripts/special-matrix/gallery.m, scripts/special-matrix/invhilb.m, scripts/statistics/base/quantile.m, scripts/statistics/tests/kruskal_wallis_test.m, scripts/statistics/tests/manova.m, scripts/statistics/tests/mcnemar_test.m: Add @nospell macro around proper names in documentation.
author Rik <rik@octave.org>
date Thu, 21 Aug 2014 15:57:18 -0700
parents 2a9114104271
children 4197fc428c7d
line wrap: on
line source

## Copyright (C) 1996-2013 Kurt Hornik
##
## This file is part of Octave.
##
## Octave 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.
##
## Octave 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 Octave; see the file COPYING.  If not, see
## <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn {Function File} {} manova (@var{x}, @var{g})
## Perform a one-way multivariate analysis of variance (MANOVA).  The
## goal is to test whether the p-dimensional population means of data
## taken from @var{k} different groups are all equal.  All data are
## assumed drawn independently from p-dimensional normal distributions
## with the same covariance matrix.
##
## The data matrix is given by @var{x}.  As usual, rows are observations
## and columns are variables.  The vector @var{g} specifies the
## corresponding group labels (e.g., numbers from 1 to @var{k}).
##
## The LR test statistic (@nospell{Wilks' Lambda}) and approximate p-values are
## computed and displayed.
## @end deftypefn

## The Hotelling-Lawley and Pillai-Bartlett test statistics are coded.
## However, they are currently disabled until they can be verified by someone
## with sufficient understanding of the algorithms.  Please feel free to
## improve this.

## Author: TF <Thomas.Fuereder@ci.tuwien.ac.at>
## Adapted-By: KH <Kurt.Hornik@wu-wien.ac.at>
## Description: One-way multivariate analysis of variance (MANOVA)

function manova (x, g)

  if (nargin != 2)
    print_usage ();
  endif

  if (isvector (x))
    error ("manova: X must not be a vector");
  endif

  [n, p] = size (x);

  if (!isvector (g) || (length (g) != n))
    error ("manova: G must be a vector of length rows (X)");
  endif

  s = sort (g);
  i = find (s (2:n) > s(1:(n-1)));
  k = length (i) + 1;

  if (k == 1)
    error ("manova: there should be at least 2 groups");
  else
    group_label = s ([1, (reshape (i, 1, k - 1) + 1)]);
  endif

  x = x - ones (n, 1) * mean (x);
  SST = x' * x;

  s = zeros (1, p);
  SSB = zeros (p, p);
  for i = 1 : k;
    v = x (find (g == group_label (i)), :);
    s = sum (v);
    SSB = SSB + s' * s / rows (v);
  endfor
  n_b = k - 1;

  SSW = SST - SSB;
  n_w = n - k;

  l = real (eig (SSB / SSW));

  if (isa (l, "single"))
    l (l < eps ("single")) = 0;
  else
    l (l < eps) = 0;
  endif

  ## Wilks' Lambda
  ## =============

  Lambda = prod (1 ./ (1 + l));

  delta = n_w + n_b - (p + n_b + 1) / 2;
  df_num = p * n_b;
  W_pval_1 = 1 - chi2cdf (- delta * log (Lambda), df_num);

  if (p < 3)
    eta = p;
  else
    eta = sqrt ((p^2 * n_b^2 - 4) / (p^2 + n_b^2 - 5));
  endif

  df_den = delta * eta - df_num / 2 + 1;

  WT = exp (- log (Lambda) / eta) - 1;
  W_pval_2 = 1 - fcdf (WT * df_den / df_num, df_num, df_den);

  if (0)

    ## Hotelling-Lawley Test
    ## =====================

    HL = sum (l);

    theta = min (p, n_b);
    u = (abs (p - n_b) - 1) / 2;
    v = (n_w - p - 1) / 2;

    df_num = theta * (2 * u + theta + 1);
    df_den = 2 * (theta * v + 1);

    HL_pval = 1 - fcdf (HL * df_den / df_num, df_num, df_den);

    ## Pillai-Bartlett
    ## ===============

    PB = sum (l ./ (1 + l));

    df_den = theta * (2 * v + theta + 1);
    PB_pval = 1 - fcdf (PB * df_den / df_num, df_num, df_den);

    printf ("\n");
    printf ("One-way MANOVA Table:\n");
    printf ("\n");
    printf ("Test             Test Statistic      Approximate p\n");
    printf ("**************************************************\n");
    printf ("Wilks            %10.4f           %10.9f \n", Lambda, W_pval_1);
    printf ("                                      %10.9f \n", W_pval_2);
    printf ("Hotelling-Lawley %10.4f           %10.9f \n", HL, HL_pval);
    printf ("Pillai-Bartlett  %10.4f           %10.9f \n", PB, PB_pval);
    printf ("\n");

  endif

  printf ("\n");
  printf ("MANOVA Results:\n");
  printf ("\n");
  printf ("# of groups:    %d\n", k);
  printf ("# of samples:   %d\n", n);
  printf ("# of variables: %d\n", p);
  printf ("\n");
  printf ("Wilks' Lambda:  %5.4f\n", Lambda);
  printf ("Approximate p:  %10.9f (chisquare approximation)\n", W_pval_1);
  printf ("                 %10.9f (F approximation)\n", W_pval_2);
  printf ("\n");

endfunction