Mercurial > octave
annotate scripts/statistics/tests/manova.m @ 21634:96518f623c91
Backed out changeset dcf8922b724b
author | Mike Miller <mtmiller@octave.org> |
---|---|
date | Wed, 20 Apr 2016 11:06:03 -0700 |
parents | dcf8922b724b |
children | bac0d6f07a3e |
rev | line source |
---|---|
19697
4197fc428c7d
maint: Update copyright notices for 2015.
John W. Eaton <jwe@octave.org>
parents:
19040
diff
changeset
|
1 ## Copyright (C) 1996-2015 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 -*- |
20852
516bb87ea72e
2015 Code Sprint: remove class of function from docstring for all m-files.
Rik <rik@octave.org>
parents:
20231
diff
changeset
|
20 ## @deftypefn {} {} manova (@var{x}, @var{g}) |
20174
d9341b422488
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
21 ## Perform a one-way multivariate analysis of variance (MANOVA). |
3200 | 22 ## |
20174
d9341b422488
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
23 ## The goal is to test whether the p-dimensional population means of data |
d9341b422488
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
24 ## taken from @var{k} different groups are all equal. All data are assumed |
d9341b422488
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
25 ## drawn independently from p-dimensional normal distributions with the same |
d9341b422488
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
26 ## covariance matrix. |
d9341b422488
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
27 ## |
d9341b422488
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
28 ## The data matrix is given by @var{x}. As usual, rows are observations and |
d9341b422488
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
29 ## columns are variables. The vector @var{g} specifies the corresponding |
d9341b422488
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
30 ## group labels (e.g., numbers from 1 to @var{k}). |
3200 | 31 ## |
19040
0850b5212619
doc: Add @nospell macro around proper names in documentation.
Rik <rik@octave.org>
parents:
17785
diff
changeset
|
32 ## The LR test statistic (@nospell{Wilks' Lambda}) and approximate p-values are |
3200 | 33 ## computed and displayed. |
20174
d9341b422488
doc: Update more docstrings to have one sentence summary as first line.
Rik <rik@octave.org>
parents:
19833
diff
changeset
|
34 ## @seealso{anova} |
3454 | 35 ## @end deftypefn |
3200 | 36 |
17785
2a9114104271
manova.m: update of description for 2 disabled tests and 2 error msg fixes
Andreas Weber <andy.weber.aw@gmail.com>
parents:
17744
diff
changeset
|
37 ## The Hotelling-Lawley and Pillai-Bartlett test statistics are coded. |
2a9114104271
manova.m: update of description for 2 disabled tests and 2 error msg fixes
Andreas Weber <andy.weber.aw@gmail.com>
parents:
17744
diff
changeset
|
38 ## However, they are currently disabled until they can be verified by someone |
2a9114104271
manova.m: update of description for 2 disabled tests and 2 error msg fixes
Andreas Weber <andy.weber.aw@gmail.com>
parents:
17744
diff
changeset
|
39 ## with sufficient understanding of the algorithms. Please feel free to |
2a9114104271
manova.m: update of description for 2 disabled tests and 2 error msg fixes
Andreas Weber <andy.weber.aw@gmail.com>
parents:
17744
diff
changeset
|
40 ## improve this. |
3426 | 41 |
3456 | 42 ## Author: TF <Thomas.Fuereder@ci.tuwien.ac.at> |
5428 | 43 ## Adapted-By: KH <Kurt.Hornik@wu-wien.ac.at> |
3456 | 44 ## Description: One-way multivariate analysis of variance (MANOVA) |
3200 | 45 |
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10680
diff
changeset
|
46 function manova (x, g) |
3200 | 47 |
48 if (nargin != 2) | |
6046 | 49 print_usage (); |
3200 | 50 endif |
51 | |
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10680
diff
changeset
|
52 if (isvector (x)) |
17785
2a9114104271
manova.m: update of description for 2 disabled tests and 2 error msg fixes
Andreas Weber <andy.weber.aw@gmail.com>
parents:
17744
diff
changeset
|
53 error ("manova: X must not be a vector"); |
3200 | 54 endif |
55 | |
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10680
diff
changeset
|
56 [n, p] = size (x); |
3200 | 57 |
19833
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
58 if (! isvector (g) || (length (g) != n)) |
17785
2a9114104271
manova.m: update of description for 2 disabled tests and 2 error msg fixes
Andreas Weber <andy.weber.aw@gmail.com>
parents:
17744
diff
changeset
|
59 error ("manova: G must be a vector of length rows (X)"); |
3200 | 60 endif |
61 | |
62 s = sort (g); | |
63 i = find (s (2:n) > s(1:(n-1))); | |
64 k = length (i) + 1; | |
3426 | 65 |
3200 | 66 if (k == 1) |
3456 | 67 error ("manova: there should be at least 2 groups"); |
3200 | 68 else |
3273 | 69 group_label = s ([1, (reshape (i, 1, k - 1) + 1)]); |
3200 | 70 endif |
71 | |
20231
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20174
diff
changeset
|
72 x -= ones (n, 1) * mean (x); |
11469
c776f063fefe
Overhaul m-script files to use common variable name between code and documentation.
Rik <octave@nomad.inbox5.com>
parents:
10680
diff
changeset
|
73 SST = x' * x; |
3200 | 74 |
75 s = zeros (1, p); | |
76 SSB = zeros (p, p); | |
77 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
|
78 v = x (find (g == group_label (i)), :); |
3200 | 79 s = sum (v); |
20231
83792dd9bcc1
Use in-place operators in m-files where possible.
Rik <rik@octave.org>
parents:
20174
diff
changeset
|
80 SSB += s' * s / rows (v); |
3200 | 81 endfor |
82 n_b = k - 1; | |
3426 | 83 |
3200 | 84 SSW = SST - SSB; |
85 n_w = n - k; | |
86 | |
87 l = real (eig (SSB / SSW)); | |
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
88 |
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
89 if (isa (l, "single")) |
19833
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
90 l(l < eps ("single")) = 0; |
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
91 else |
19833
9fc020886ae9
maint: Clean up m-files to follow Octave coding conventions.
Rik <rik@octave.org>
parents:
19697
diff
changeset
|
92 l(l < eps) = 0; |
7795
df9519e9990c
Handle single precision eps values
David Bateman <dbateman@free.fr>
parents:
7017
diff
changeset
|
93 endif |
3200 | 94 |
95 ## Wilks' Lambda | |
96 ## ============= | |
97 | |
98 Lambda = prod (1 ./ (1 + l)); | |
3426 | 99 |
11589
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
100 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
|
101 df_num = p * n_b; |
10680
e00de2d5263c
Replace calls to obsolete chisquare_cdf with chi2cdf.
Rik <octave@nomad.inbox5.com>
parents:
9245
diff
changeset
|
102 W_pval_1 = 1 - chi2cdf (- delta * log (Lambda), df_num); |
3426 | 103 |
3200 | 104 if (p < 3) |
105 eta = p; | |
106 else | |
11589
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
107 eta = sqrt ((p^2 * n_b^2 - 4) / (p^2 + n_b^2 - 5)); |
3200 | 108 endif |
109 | |
11589
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
110 df_den = delta * eta - df_num / 2 + 1; |
3426 | 111 |
11589
b0084095098e
missing semicolons in script files
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
112 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
|
113 W_pval_2 = 1 - fcdf (WT * df_den / df_num, df_num, df_den); |
3200 | 114 |
115 if (0) | |
116 | |
117 ## Hotelling-Lawley Test | |
118 ## ===================== | |
3426 | 119 |
3200 | 120 HL = sum (l); |
3426 | 121 |
3200 | 122 theta = min (p, n_b); |
3426 | 123 u = (abs (p - n_b) - 1) / 2; |
3200 | 124 v = (n_w - p - 1) / 2; |
125 | |
126 df_num = theta * (2 * u + theta + 1); | |
127 df_den = 2 * (theta * v + 1); | |
128 | |
13752
6f068e3f3f9c
Change f_cdf references to fcdf in statistics/test directory (Bug #34628)
Rik <octave@nomad.inbox5.com>
parents:
11589
diff
changeset
|
129 HL_pval = 1 - fcdf (HL * df_den / df_num, df_num, df_den); |
3200 | 130 |
131 ## Pillai-Bartlett | |
132 ## =============== | |
3426 | 133 |
3200 | 134 PB = sum (l ./ (1 + l)); |
135 | |
136 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
|
137 PB_pval = 1 - fcdf (PB * df_den / df_num, df_num, df_den); |
3200 | 138 |
21634
96518f623c91
Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents:
21633
diff
changeset
|
139 printf ("\n"); |
96518f623c91
Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents:
21633
diff
changeset
|
140 printf ("One-way MANOVA Table:\n"); |
96518f623c91
Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents:
21633
diff
changeset
|
141 printf ("\n"); |
96518f623c91
Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents:
21633
diff
changeset
|
142 printf ("Test Test Statistic Approximate p\n"); |
96518f623c91
Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents:
21633
diff
changeset
|
143 printf ("**************************************************\n"); |
96518f623c91
Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents:
21633
diff
changeset
|
144 printf ("Wilks %10.4f %10.9f \n", Lambda, W_pval_1); |
96518f623c91
Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents:
21633
diff
changeset
|
145 printf (" %10.9f \n", W_pval_2); |
96518f623c91
Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents:
21633
diff
changeset
|
146 printf ("Hotelling-Lawley %10.4f %10.9f \n", HL, HL_pval); |
96518f623c91
Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents:
21633
diff
changeset
|
147 printf ("Pillai-Bartlett %10.4f %10.9f \n", PB, PB_pval); |
96518f623c91
Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents:
21633
diff
changeset
|
148 printf ("\n"); |
3200 | 149 |
150 endif | |
151 | |
21634
96518f623c91
Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents:
21633
diff
changeset
|
152 printf ("\n"); |
96518f623c91
Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents:
21633
diff
changeset
|
153 printf ("MANOVA Results:\n"); |
96518f623c91
Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents:
21633
diff
changeset
|
154 printf ("\n"); |
96518f623c91
Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents:
21633
diff
changeset
|
155 printf ("# of groups: %d\n", k); |
96518f623c91
Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents:
21633
diff
changeset
|
156 printf ("# of samples: %d\n", n); |
96518f623c91
Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents:
21633
diff
changeset
|
157 printf ("# of variables: %d\n", p); |
96518f623c91
Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents:
21633
diff
changeset
|
158 printf ("\n"); |
96518f623c91
Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents:
21633
diff
changeset
|
159 printf ("Wilks' Lambda: %5.4f\n", Lambda); |
96518f623c91
Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents:
21633
diff
changeset
|
160 printf ("Approximate p: %10.9f (chisquare approximation)\n", W_pval_1); |
96518f623c91
Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents:
21633
diff
changeset
|
161 printf (" %10.9f (F approximation)\n", W_pval_2); |
96518f623c91
Backed out changeset dcf8922b724b
Mike Miller <mtmiller@octave.org>
parents:
21633
diff
changeset
|
162 printf ("\n"); |
3426 | 163 |
3200 | 164 endfunction |
17338
1c89599167a6
maint: End m-files with 1 blank line.
Rik <rik@octave.org>
parents:
14138
diff
changeset
|
165 |