Mercurial > octave-nkf
annotate scripts/statistics/tests/manova.m @ 14138:72c96de7a403 stable
maint: update copyright notices for 2012
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Mon, 02 Jan 2012 14:25:41 -0500 |
| parents | 6f068e3f3f9c |
| children | 1c89599167a6 |
| rev | line source |
|---|---|
|
14138
72c96de7a403
maint: update copyright notices for 2012
John W. Eaton <jwe@octave.org>
parents:
13752
diff
changeset
|
1 ## Copyright (C) 1996-2012 Kurt Hornik |
| 3200 | 2 ## |
| 3922 | 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 | |
| 7016 | 7 ## the Free Software Foundation; either version 3 of the License, or (at |
| 8 ## your option) any later version. | |
| 3200 | 9 ## |
| 3922 | 10 ## Octave is distributed in the hope that it will be useful, but |
| 3200 | 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 | |
| 7016 | 16 ## along with Octave; see the file COPYING. If not, see |
| 17 ## <http://www.gnu.org/licenses/>. | |
| 3200 | 18 |
| 3454 | 19 ## -*- texinfo -*- |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10680
diff
changeset
|
20 ## @deftypefn {Function File} {} manova (@var{x}, @var{g}) |
|
9051
1bf0ce0930be
Grammar check TexInfo in all .m files
Rik <rdrider0-list@yahoo.com>
parents:
8920
diff
changeset
|
21 ## Perform a one-way multivariate analysis of variance (MANOVA). The |
| 3200 | 22 ## goal is to test whether the p-dimensional population means of data |
| 3454 | 23 ## taken from @var{k} different groups are all equal. All data are |
| 24 ## assumed drawn independently from p-dimensional normal distributions | |
| 25 ## with the same covariance matrix. | |
| 3200 | 26 ## |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10680
diff
changeset
|
27 ## The data matrix is given by @var{x}. As usual, rows are observations |
| 3454 | 28 ## and columns are variables. The vector @var{g} specifies the |
| 29 ## corresponding group labels (e.g., numbers from 1 to @var{k}). | |
| 3200 | 30 ## |
| 31 ## The LR test statistic (Wilks' Lambda) and approximate p-values are | |
| 32 ## computed and displayed. | |
| 3454 | 33 ## @end deftypefn |
| 3200 | 34 |
| 35 ## Three test statistics (Wilks, Hotelling-Lawley, and Pillai-Bartlett) | |
| 36 ## and corresponding approximate p-values are calculated and displayed. | |
|
13752
6f068e3f3f9c
Change f_cdf references to fcdf in statistics/test directory (Bug #34628)
Rik <octave@nomad.inbox5.com>
parents:
11589
diff
changeset
|
37 ## (Currently NOT because the fcdf respectively betai code is too bad.) |
| 3426 | 38 |
| 3456 | 39 ## Author: TF <Thomas.Fuereder@ci.tuwien.ac.at> |
| 5428 | 40 ## Adapted-By: KH <Kurt.Hornik@wu-wien.ac.at> |
| 3456 | 41 ## Description: One-way multivariate analysis of variance (MANOVA) |
| 3200 | 42 |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10680
diff
changeset
|
43 function manova (x, g) |
| 3200 | 44 |
| 45 if (nargin != 2) | |
| 6046 | 46 print_usage (); |
| 3200 | 47 endif |
| 48 | |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10680
diff
changeset
|
49 if (isvector (x)) |
| 3456 | 50 error ("manova: Y must not be a vector"); |
| 3200 | 51 endif |
| 52 | |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10680
diff
changeset
|
53 [n, p] = size (x); |
| 3200 | 54 |
| 4030 | 55 if (!isvector (g) || (length (g) != n)) |
|
11472
1740012184f9
Use uppercase for variable names in error() strings to match Info documentation. Only m-files done.
Rik <octave@nomad.inbox5.com>
parents:
11469
diff
changeset
|
56 error ("manova: G must be a vector of length rows (Y)"); |
| 3200 | 57 endif |
| 58 | |
| 59 s = sort (g); | |
| 60 i = find (s (2:n) > s(1:(n-1))); | |
| 61 k = length (i) + 1; | |
| 3426 | 62 |
| 3200 | 63 if (k == 1) |
| 3456 | 64 error ("manova: there should be at least 2 groups"); |
| 3200 | 65 else |
| 3273 | 66 group_label = s ([1, (reshape (i, 1, k - 1) + 1)]); |
| 3200 | 67 endif |
| 68 | |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10680
diff
changeset
|
69 x = x - ones (n, 1) * mean (x); |
|
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10680
diff
changeset
|
70 SST = x' * x; |
| 3200 | 71 |
| 72 s = zeros (1, p); | |
| 73 SSB = zeros (p, p); | |
| 74 for i = 1 : k; | |
|
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10680
diff
changeset
|
75 v = x (find (g == group_label (i)), :); |
| 3200 | 76 s = sum (v); |
| 77 SSB = SSB + s' * s / rows (v); | |
| 78 endfor | |
| 79 n_b = k - 1; | |
| 3426 | 80 |
| 3200 | 81 SSW = SST - SSB; |
| 82 n_w = n - k; | |
| 83 | |
| 84 l = real (eig (SSB / SSW)); | |
|
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
85 |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
86 if (isa (l, "single")) |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
87 l (l < eps ("single")) = 0; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
88 else |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
89 l (l < eps) = 0; |
|
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
90 endif |
| 3200 | 91 |
| 92 ## Wilks' Lambda | |
| 93 ## ============= | |
| 94 | |
| 95 Lambda = prod (1 ./ (1 + l)); | |
| 3426 | 96 |
|
11589
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
97 delta = n_w + n_b - (p + n_b + 1) / 2; |
|
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
98 df_num = p * n_b; |
|
10680
e00de2d5263c
Replace calls to obsolete chisquare_cdf with chi2cdf.
Rik <octave@nomad.inbox5.com>
parents:
9245
diff
changeset
|
99 W_pval_1 = 1 - chi2cdf (- delta * log (Lambda), df_num); |
| 3426 | 100 |
| 3200 | 101 if (p < 3) |
| 102 eta = p; | |
| 103 else | |
|
11589
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
104 eta = sqrt ((p^2 * n_b^2 - 4) / (p^2 + n_b^2 - 5)); |
| 3200 | 105 endif |
| 106 | |
|
11589
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
107 df_den = delta * eta - df_num / 2 + 1; |
| 3426 | 108 |
|
11589
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
109 WT = exp (- log (Lambda) / eta) - 1; |
|
13752
6f068e3f3f9c
Change f_cdf references to fcdf in statistics/test directory (Bug #34628)
Rik <octave@nomad.inbox5.com>
parents:
11589
diff
changeset
|
110 W_pval_2 = 1 - fcdf (WT * df_den / df_num, df_num, df_den); |
| 3200 | 111 |
| 112 if (0) | |
| 113 | |
| 114 ## Hotelling-Lawley Test | |
| 115 ## ===================== | |
| 3426 | 116 |
| 3200 | 117 HL = sum (l); |
| 3426 | 118 |
| 3200 | 119 theta = min (p, n_b); |
| 3426 | 120 u = (abs (p - n_b) - 1) / 2; |
| 3200 | 121 v = (n_w - p - 1) / 2; |
| 122 | |
| 123 df_num = theta * (2 * u + theta + 1); | |
| 124 df_den = 2 * (theta * v + 1); | |
| 125 | |
|
13752
6f068e3f3f9c
Change f_cdf references to fcdf in statistics/test directory (Bug #34628)
Rik <octave@nomad.inbox5.com>
parents:
11589
diff
changeset
|
126 HL_pval = 1 - fcdf (HL * df_den / df_num, df_num, df_den); |
| 3200 | 127 |
| 128 ## Pillai-Bartlett | |
| 129 ## =============== | |
| 3426 | 130 |
| 3200 | 131 PB = sum (l ./ (1 + l)); |
| 132 | |
| 133 df_den = theta * (2 * v + theta + 1); | |
|
13752
6f068e3f3f9c
Change f_cdf references to fcdf in statistics/test directory (Bug #34628)
Rik <octave@nomad.inbox5.com>
parents:
11589
diff
changeset
|
134 PB_pval = 1 - fcdf (PB * df_den / df_num, df_num, df_den); |
| 3200 | 135 |
| 136 printf ("\n"); | |
| 137 printf ("One-way MANOVA Table:\n"); | |
| 3426 | 138 printf ("\n"); |
| 3200 | 139 printf ("Test Test Statistic Approximate p\n"); |
| 140 printf ("**************************************************\n"); | |
| 141 printf ("Wilks %10.4f %10.9f \n", Lambda, W_pval_1); | |
| 142 printf (" %10.9f \n", W_pval_2); | |
| 143 printf ("Hotelling-Lawley %10.4f %10.9f \n", HL, HL_pval); | |
| 144 printf ("Pillai-Bartlett %10.4f %10.9f \n", PB, PB_pval); | |
| 145 printf ("\n"); | |
| 146 | |
| 147 endif | |
| 148 | |
| 149 printf ("\n"); | |
| 150 printf ("MANOVA Results:\n"); | |
| 151 printf ("\n"); | |
| 3456 | 152 printf ("# of groups: %d\n", k); |
| 153 printf ("# of samples: %d\n", n); | |
| 154 printf ("# of variables: %d\n", p); | |
| 3426 | 155 printf ("\n"); |
| 3456 | 156 printf ("Wilks' Lambda: %5.4f\n", Lambda); |
| 157 printf ("Approximate p: %10.9f (chisquare approximation)\n", W_pval_1); | |
| 3200 | 158 printf (" %10.9f (F approximation)\n", W_pval_2); |
| 159 printf ("\n"); | |
| 3426 | 160 |
| 3200 | 161 endfunction |