annotate scripts/statistics/tests/cor_test.m @ 7017:a1dbe9d80eee

[project @ 2007-10-12 21:27:11 by jwe]
author jwe
date Fri, 12 Oct 2007 21:27:37 +0000
parents 93c65f2a5668
children e07e93c04080
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
7017
a1dbe9d80eee [project @ 2007-10-12 21:27:11 by jwe]
jwe
parents: 7016
diff changeset
1 ## Copyright (C) 1995, 1996, 1997, 1998, 2000, 2002, 2005, 2006, 2007
a1dbe9d80eee [project @ 2007-10-12 21:27:11 by jwe]
jwe
parents: 7016
diff changeset
2 ## Kurt Hornik
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3200
diff changeset
3 ##
3922
38c61cbf086c [project @ 2002-05-01 06:48:35 by jwe]
jwe
parents: 3457
diff changeset
4 ## This file is part of Octave.
38c61cbf086c [project @ 2002-05-01 06:48:35 by jwe]
jwe
parents: 3457
diff changeset
5 ##
38c61cbf086c [project @ 2002-05-01 06:48:35 by jwe]
jwe
parents: 3457
diff changeset
6 ## Octave is free software; you can redistribute it and/or modify it
38c61cbf086c [project @ 2002-05-01 06:48:35 by jwe]
jwe
parents: 3457
diff changeset
7 ## under the terms of the GNU General Public License as published by
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6046
diff changeset
8 ## the Free Software Foundation; either version 3 of the License, or (at
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6046
diff changeset
9 ## your option) any later version.
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3200
diff changeset
10 ##
3922
38c61cbf086c [project @ 2002-05-01 06:48:35 by jwe]
jwe
parents: 3457
diff changeset
11 ## Octave is distributed in the hope that it will be useful, but
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3200
diff changeset
14 ## General Public License for more details.
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3200
diff changeset
15 ##
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
16 ## You should have received a copy of the GNU General Public License
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6046
diff changeset
17 ## along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6046
diff changeset
18 ## <http://www.gnu.org/licenses/>.
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
19
3454
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
20 ## -*- texinfo -*-
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
21 ## @deftypefn {Function File} {} cor_test (@var{x}, @var{y}, @var{alt}, @var{method})
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
22 ## Test whether two samples @var{x} and @var{y} come from uncorrelated
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
23 ## populations.
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
24 ##
3454
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
25 ## The optional argument string @var{alt} describes the alternative
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
26 ## hypothesis, and can be @code{"!="} or @code{"<>"} (non-zero),
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
27 ## @code{">"} (greater than 0), or @code{"<"} (less than 0). The
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
28 ## default is the two-sided case.
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
29 ##
3454
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
30 ## The optional argument string @var{method} specifies on which
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
31 ## correlation coefficient the test should be based. If @var{method} is
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
32 ## @code{"pearson"} (default), the (usual) Pearson's product moment
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
33 ## correlation coefficient is used. In this case, the data should come
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
34 ## from a bivariate normal distribution. Otherwise, the other two
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
35 ## methods offer nonparametric alternatives. If @var{method} is
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
36 ## @code{"kendall"}, then Kendall's rank correlation tau is used. If
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
37 ## @var{method} is @code{"spearman"}, then Spearman's rank correlation
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
38 ## rho is used. Only the first character is necessary.
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
39 ##
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
40 ## The output is a structure with the following elements:
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
41 ##
3454
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
42 ## @table @var
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
43 ## @item pval
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
44 ## The p-value of the test.
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
45 ## @item stat
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
46 ## The value of the test statistic.
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
47 ## @item dist
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
48 ## The distribution of the test statistic.
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
49 ## @item params
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
50 ## The parameters of the null distribution of the test statistic.
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
51 ## @item alternative
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
52 ## The alternative hypothesis.
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
53 ## @item method
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
54 ## The method used for testing.
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
55 ## @end table
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
56 ##
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
57 ## If no output argument is given, the p-value is displayed.
d8b731d3f7a3 [project @ 2000-01-18 10:13:31 by jwe]
jwe
parents: 3426
diff changeset
58 ## @end deftypefn
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
59
3456
434790acb067 [project @ 2000-01-19 06:58:51 by jwe]
jwe
parents: 3454
diff changeset
60 ## Author: FL <Friedrich.Leisch@ci.tuwien.ac.at>
5428
2a16423e4aa0 [project @ 2005-08-23 18:38:27 by jwe]
jwe
parents: 5307
diff changeset
61 ## Adapted-by: KH <Kurt.Hornik@wu-wien.ac.at>
3456
434790acb067 [project @ 2000-01-19 06:58:51 by jwe]
jwe
parents: 3454
diff changeset
62 ## Description: Test for zero correlation
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
63
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
64 function t = cor_test (X, Y, ALTERNATIVE, METHOD)
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3200
diff changeset
65
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
66 if ((nargin < 2) || (nargin > 4))
6046
34f96dd5441b [project @ 2006-10-10 16:10:25 by jwe]
jwe
parents: 5443
diff changeset
67 print_usage ();
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
68 endif
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
69
4030
22bd65326ec1 [project @ 2002-08-09 18:58:13 by jwe]
jwe
parents: 3922
diff changeset
70 if (!isvector (X) || !isvector (Y) || length (X) != length (Y))
3456
434790acb067 [project @ 2000-01-19 06:58:51 by jwe]
jwe
parents: 3454
diff changeset
71 error ("cor_test: X and Y must be vectors of the same length")
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
72 endif
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
73
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
74 if (nargin < 3)
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
75 ALTERNATIVE = "!=";
5443
ec8c33dcd1bf [project @ 2005-09-08 01:40:57 by jwe]
jwe
parents: 5428
diff changeset
76 elseif (! ischar (ALTERNATIVE))
3456
434790acb067 [project @ 2000-01-19 06:58:51 by jwe]
jwe
parents: 3454
diff changeset
77 error ("cor_test: ALTERNATIVE must be a string");
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
78 endif
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
79
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
80 if (nargin < 4)
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
81 METHOD = "pearson";
5443
ec8c33dcd1bf [project @ 2005-09-08 01:40:57 by jwe]
jwe
parents: 5428
diff changeset
82 elseif (! ischar (METHOD))
3456
434790acb067 [project @ 2000-01-19 06:58:51 by jwe]
jwe
parents: 3454
diff changeset
83 error ("cor_test: METHOD must be a string");
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
84 endif
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
85
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
86 n = length (X);
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
87 m = METHOD (1);
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
88
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
89 if (m == "p")
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
90 r = cor (X, Y);
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
91 df = n - 2;
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
92 t.method = "Pearson's product moment correlation";
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
93 t.params = df;
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
94 t.stat = sqrt (df) .* r / sqrt (1 - r.^2);
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
95 t.dist = "t";
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
96 cdf = t_cdf (t.stat, df);
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
97 elseif (m == "k")
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
98 tau = kendall (X, Y);
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3200
diff changeset
99 t.method = "Kendall's rank correlation tau";
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
100 t.params = [];
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
101 t.stat = tau / sqrt ((2 * (2*n+5)) / (9*n*(n-1)));
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
102 t.dist = "stdnormal";
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
103 cdf = stdnormal_cdf (t.stat);
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
104 elseif (m == "s")
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
105 rho = spearman (X, Y);
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
106 t.method = "Spearman's rank correlation rho";
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
107 t.params = [];
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
108 t.stat = sqrt (n-1) * (rho - 6/(n^3-n));
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3200
diff changeset
109 t.dist = "stdnormal";
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
110 cdf = stdnormal_cdf (t.stat);
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
111 else
3456
434790acb067 [project @ 2000-01-19 06:58:51 by jwe]
jwe
parents: 3454
diff changeset
112 error ("cor_test: method `%s' not recognized", METHOD)
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
113 endif
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
114
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
115 if (strcmp (ALTERNATIVE, "!=") || strcmp (ALTERNATIVE, "<>"))
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
116 t.pval = 2 * min (cdf, 1 - cdf);
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
117 elseif (strcmp (ALTERNATIVE, ">"))
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
118 t.pval = 1 - cdf;
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
119 elseif (strcmp (ALTERNATIVE, "<"))
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
120 t.pval = cdf;
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
121 else
3456
434790acb067 [project @ 2000-01-19 06:58:51 by jwe]
jwe
parents: 3454
diff changeset
122 error ("cor_test: alternative `%s' not recognized", ALTERNATIVE);
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
123 endif
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
124
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
125 t.alternative = ALTERNATIVE;
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
126
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
127 if (nargout == 0)
3456
434790acb067 [project @ 2000-01-19 06:58:51 by jwe]
jwe
parents: 3454
diff changeset
128 printf ("pval: %g\n", t.pval);
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
129 endif
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3200
diff changeset
130
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
131 endfunction