Mercurial > hg > camir-aes2014

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

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