Mercurial > octave-nkf
annotate scripts/statistics/tests/manova.m @ 7795:df9519e9990c
Handle single precision eps values
| author | David Bateman <dbateman@free.fr> |
|---|---|
| date | Mon, 12 May 2008 22:57:11 +0200 |
| parents | a1dbe9d80eee |
| children | eb63fbe60fab |
| rev | line source |
|---|---|
| 7017 | 1 ## Copyright (C) 1996, 1997, 1998, 1999, 2000, 2002, 2005, 2006, 2007 |
| 2 ## Kurt Hornik | |
| 3200 | 3 ## |
| 3922 | 4 ## This file is part of Octave. |
| 5 ## | |
| 6 ## Octave is free software; you can redistribute it and/or modify it | |
| 7 ## under the terms of the GNU General Public License as published by | |
| 7016 | 8 ## the Free Software Foundation; either version 3 of the License, or (at |
| 9 ## your option) any later version. | |
| 3200 | 10 ## |
| 3922 | 11 ## Octave is distributed in the hope that it will be useful, but |
| 3200 | 12 ## WITHOUT ANY WARRANTY; without even the implied warranty of |
| 13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |
| 14 ## General Public License for more details. | |
| 15 ## | |
| 16 ## You should have received a copy of the GNU General Public License | |
| 7016 | 17 ## along with Octave; see the file COPYING. If not, see |
| 18 ## <http://www.gnu.org/licenses/>. | |
| 3200 | 19 |
| 3454 | 20 ## -*- texinfo -*- |
| 21 ## @deftypefn {Function File} {} manova (@var{y}, @var{g}) | |
| 22 ## Perform a one-way multivariate analysis of variance (MANOVA). The | |
| 3200 | 23 ## goal is to test whether the p-dimensional population means of data |
| 3454 | 24 ## taken from @var{k} different groups are all equal. All data are |
| 25 ## assumed drawn independently from p-dimensional normal distributions | |
| 26 ## with the same covariance matrix. | |
| 3200 | 27 ## |
| 3454 | 28 ## The data matrix is given by @var{y}. As usual, rows are observations |
| 29 ## and columns are variables. The vector @var{g} specifies the | |
| 30 ## corresponding group labels (e.g., numbers from 1 to @var{k}). | |
| 3200 | 31 ## |
| 32 ## The LR test statistic (Wilks' Lambda) and approximate p-values are | |
| 33 ## computed and displayed. | |
| 3454 | 34 ## @end deftypefn |
| 3200 | 35 |
| 36 ## Three test statistics (Wilks, Hotelling-Lawley, and Pillai-Bartlett) | |
| 37 ## and corresponding approximate p-values are calculated and displayed. | |
| 38 ## (Currently NOT because the f_cdf respectively betai code is too bad.) | |
| 3426 | 39 |
| 3456 | 40 ## Author: TF <Thomas.Fuereder@ci.tuwien.ac.at> |
| 5428 | 41 ## Adapted-By: KH <Kurt.Hornik@wu-wien.ac.at> |
| 3456 | 42 ## Description: One-way multivariate analysis of variance (MANOVA) |
| 3200 | 43 |
| 44 function manova (Y, g) | |
| 45 | |
| 46 if (nargin != 2) | |
| 6046 | 47 print_usage (); |
| 3200 | 48 endif |
| 49 | |
| 4030 | 50 if (isvector (Y)) |
| 3456 | 51 error ("manova: Y must not be a vector"); |
| 3200 | 52 endif |
| 53 | |
| 54 [n, p] = size (Y); | |
| 55 | |
| 4030 | 56 if (!isvector (g) || (length (g) != n)) |
| 3456 | 57 error ("manova: g must be a vector of length rows (Y)"); |
| 3200 | 58 endif |
| 59 | |
| 60 s = sort (g); | |
| 61 i = find (s (2:n) > s(1:(n-1))); | |
| 62 k = length (i) + 1; | |
| 3426 | 63 |
| 3200 | 64 if (k == 1) |
| 3456 | 65 error ("manova: there should be at least 2 groups"); |
| 3200 | 66 else |
| 3273 | 67 group_label = s ([1, (reshape (i, 1, k - 1) + 1)]); |
| 3200 | 68 endif |
| 69 | |
| 70 Y = Y - ones (n, 1) * mean (Y); | |
| 71 SST = Y' * Y; | |
| 72 | |
| 73 s = zeros (1, p); | |
| 74 SSB = zeros (p, p); | |
| 75 for i = 1 : k; | |
| 76 v = Y (find (g == group_label (i)), :); | |
| 77 s = sum (v); | |
| 78 SSB = SSB + s' * s / rows (v); | |
| 79 endfor | |
| 80 n_b = k - 1; | |
| 3426 | 81 |
| 3200 | 82 SSW = SST - SSB; |
| 83 n_w = n - k; | |
| 84 | |
| 85 l = real (eig (SSB / SSW)); | |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
86 |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
87 if (isa (l, "single")) |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
88 l (l < eps ("single")) = 0; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
89 else |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
90 l (l < eps) = 0; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
91 endif |
| 3200 | 92 |
| 93 ## Wilks' Lambda | |
| 94 ## ============= | |
| 95 | |
| 96 Lambda = prod (1 ./ (1 + l)); | |
| 3426 | 97 |
| 3200 | 98 delta = n_w + n_b - (p + n_b + 1) / 2 |
| 99 df_num = p * n_b | |
| 100 W_pval_1 = 1 - chisquare_cdf (- delta * log (Lambda), df_num); | |
| 3426 | 101 |
| 3200 | 102 if (p < 3) |
| 103 eta = p; | |
| 104 else | |
| 105 eta = sqrt ((p^2 * n_b^2 - 4) / (p^2 + n_b^2 - 5)) | |
| 106 endif | |
| 107 | |
| 108 df_den = delta * eta - df_num / 2 + 1 | |
| 3426 | 109 |
| 3200 | 110 WT = exp (- log (Lambda) / eta) - 1 |
| 111 W_pval_2 = 1 - f_cdf (WT * df_den / df_num, df_num, df_den); | |
| 112 | |
| 113 if (0) | |
| 114 | |
| 115 ## Hotelling-Lawley Test | |
| 116 ## ===================== | |
| 3426 | 117 |
| 3200 | 118 HL = sum (l); |
| 3426 | 119 |
| 3200 | 120 theta = min (p, n_b); |
| 3426 | 121 u = (abs (p - n_b) - 1) / 2; |
| 3200 | 122 v = (n_w - p - 1) / 2; |
| 123 | |
| 124 df_num = theta * (2 * u + theta + 1); | |
| 125 df_den = 2 * (theta * v + 1); | |
| 126 | |
| 127 HL_pval = 1 - f_cdf (HL * df_den / df_num, df_num, df_den); | |
| 128 | |
| 129 ## Pillai-Bartlett | |
| 130 ## =============== | |
| 3426 | 131 |
| 3200 | 132 PB = sum (l ./ (1 + l)); |
| 133 | |
| 134 df_den = theta * (2 * v + theta + 1); | |
| 135 PB_pval = 1 - f_cdf (PB * df_den / df_num, df_num, df_den); | |
| 136 | |
| 137 printf ("\n"); | |
| 138 printf ("One-way MANOVA Table:\n"); | |
| 3426 | 139 printf ("\n"); |
| 3200 | 140 printf ("Test Test Statistic Approximate p\n"); |
| 141 printf ("**************************************************\n"); | |
| 142 printf ("Wilks %10.4f %10.9f \n", Lambda, W_pval_1); | |
| 143 printf (" %10.9f \n", W_pval_2); | |
| 144 printf ("Hotelling-Lawley %10.4f %10.9f \n", HL, HL_pval); | |
| 145 printf ("Pillai-Bartlett %10.4f %10.9f \n", PB, PB_pval); | |
| 146 printf ("\n"); | |
| 147 | |
| 148 endif | |
| 149 | |
| 150 printf ("\n"); | |
| 151 printf ("MANOVA Results:\n"); | |
| 152 printf ("\n"); | |
| 3456 | 153 printf ("# of groups: %d\n", k); |
| 154 printf ("# of samples: %d\n", n); | |
| 155 printf ("# of variables: %d\n", p); | |
| 3426 | 156 printf ("\n"); |
| 3456 | 157 printf ("Wilks' Lambda: %5.4f\n", Lambda); |
| 158 printf ("Approximate p: %10.9f (chisquare approximation)\n", W_pval_1); | |
| 3200 | 159 printf (" %10.9f (F approximation)\n", W_pval_2); |
| 160 printf ("\n"); | |
| 3426 | 161 |
| 3200 | 162 endfunction |