comparison toolboxes/FullBNT-1.0.7/KPMstats/cond_indep_fisher_z.m @ 0:e9a9cd732c1e tip

first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:e9a9cd732c1e
1 function [CI, r, p] = cond_indep_fisher_z(X, Y, S, C, N, alpha)
2 % COND_INDEP_FISHER_Z Test if X indep Y given Z using Fisher's Z test
3 % CI = cond_indep_fisher_z(X, Y, S, C, N, alpha)
4 %
5 % C is the covariance (or correlation) matrix
6 % N is the sample size
7 % alpha is the significance level (default: 0.05)
8 %
9 % See p133 of T. Anderson, "An Intro. to Multivariate Statistical Analysis", 1984
10
11 if nargin < 6, alpha = 0.05; end
12
13 r = partial_corr_coef(C, X, Y, S);
14 z = 0.5*log( (1+r)/(1-r) );
15 z0 = 0;
16 W = sqrt(N - length(S) - 3)*(z-z0); % W ~ N(0,1)
17 cutoff = norminv(1 - 0.5*alpha); % P(|W| <= cutoff) = 0.95
18 %cutoff = mynorminv(1 - 0.5*alpha); % P(|W| <= cutoff) = 0.95
19 if abs(W) < cutoff
20 CI = 1;
21 else % reject the null hypothesis that rho = 0
22 CI = 0;
23 end
24 p = normcdf(W);
25 %p = mynormcdf(W);
26
27 %%%%%%%%%
28
29 function p = normcdf(x,mu,sigma)
30 %NORMCDF Normal cumulative distribution function (cdf).
31 % P = NORMCDF(X,MU,SIGMA) computes the normal cdf with mean MU and
32 % standard deviation SIGMA at the values in X.
33 %
34 % The size of P is the common size of X, MU and SIGMA. A scalar input
35 % functions as a constant matrix of the same size as the other inputs.
36 %
37 % Default values for MU and SIGMA are 0 and 1 respectively.
38
39 % References:
40 % [1] M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
41 % Functions", Government Printing Office, 1964, 26.2.
42
43 % Copyright (c) 1993-98 by The MathWorks, Inc.
44 % $Revision: 1.1.1.1 $ $Date: 2005/04/26 02:29:18 $
45
46 if nargin < 3,
47 sigma = 1;
48 end
49
50 if nargin < 2;
51 mu = 0;
52 end
53
54 [errorcode x mu sigma] = distchck(3,x,mu,sigma);
55
56 if errorcode > 0
57 error('Requires non-scalar arguments to match in size.');
58 end
59
60 % Initialize P to zero.
61 p = zeros(size(x));
62
63 % Return NaN if SIGMA is not positive.
64 k1 = find(sigma <= 0);
65 if any(k1)
66 tmp = NaN;
67 p(k1) = tmp(ones(size(k1)));
68 end
69
70 % Express normal CDF in terms of the error function.
71 k = find(sigma > 0);
72 if any(k)
73 p(k) = 0.5 * erfc( - (x(k) - mu(k)) ./ (sigma(k) * sqrt(2)));
74 end
75
76 % Make sure that round-off errors never make P greater than 1.
77 k2 = find(p > 1);
78 if any(k2)
79 p(k2) = ones(size(k2));
80 end
81
82 %%%%%%%%
83
84 function x = norminv(p,mu,sigma);
85 %NORMINV Inverse of the normal cumulative distribution function (cdf).
86 % X = NORMINV(P,MU,SIGMA) finds the inverse of the normal cdf with
87 % mean, MU, and standard deviation, SIGMA.
88 %
89 % The size of X is the common size of the input arguments. A scalar input
90 % functions as a constant matrix of the same size as the other inputs.
91 %
92 % Default values for MU and SIGMA are 0 and 1 respectively.
93
94 % References:
95 % [1] M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
96 % Functions", Government Printing Office, 1964, 7.1.1 and 26.2.2
97
98 % Copyright (c) 1993-98 by The MathWorks, Inc.
99 % $Revision: 1.1.1.1 $ $Date: 2005/04/26 02:29:18 $
100
101 if nargin < 3,
102 sigma = 1;
103 end
104
105 if nargin < 2;
106 mu = 0;
107 end
108
109 [errorcode p mu sigma] = distchck(3,p,mu,sigma);
110
111 if errorcode > 0
112 error('Requires non-scalar arguments to match in size.');
113 end
114
115 % Allocate space for x.
116 x = zeros(size(p));
117
118 % Return NaN if the arguments are outside their respective limits.
119 k = find(sigma <= 0 | p < 0 | p > 1);
120 if any(k)
121 tmp = NaN;
122 x(k) = tmp(ones(size(k)));
123 end
124
125 % Put in the correct values when P is either 0 or 1.
126 k = find(p == 0);
127 if any(k)
128 tmp = Inf;
129 x(k) = -tmp(ones(size(k)));
130 end
131
132 k = find(p == 1);
133 if any(k)
134 tmp = Inf;
135 x(k) = tmp(ones(size(k)));
136 end
137
138 % Compute the inverse function for the intermediate values.
139 k = find(p > 0 & p < 1 & sigma > 0);
140 if any(k),
141 x(k) = sqrt(2) * sigma(k) .* erfinv(2 * p(k) - 1) + mu(k);
142 end