wolffd@0: function [CI, r, p] = cond_indep_fisher_z(X, Y, S, C, N, alpha)
wolffd@0: % COND_INDEP_FISHER_Z Test if X indep Y given Z using Fisher's Z test
wolffd@0: % CI = cond_indep_fisher_z(X, Y, S, C, N, alpha)
wolffd@0: %
wolffd@0: % C is the covariance (or correlation) matrix
wolffd@0: % N is the sample size
wolffd@0: % alpha is the significance level (default: 0.05)
wolffd@0: %
wolffd@0: % See p133 of T. Anderson, "An Intro. to Multivariate Statistical Analysis", 1984
wolffd@0: 
wolffd@0: if nargin < 6, alpha = 0.05; end
wolffd@0: 
wolffd@0: r = partial_corr_coef(C, X, Y, S);
wolffd@0: z = 0.5*log( (1+r)/(1-r) );
wolffd@0: z0 = 0;
wolffd@0: W = sqrt(N - length(S) - 3)*(z-z0); % W ~ N(0,1)
wolffd@0: cutoff = norminv(1 - 0.5*alpha); % P(|W| <= cutoff) = 0.95
wolffd@0: %cutoff = mynorminv(1 - 0.5*alpha); % P(|W| <= cutoff) = 0.95
wolffd@0: if abs(W) < cutoff
wolffd@0:   CI = 1;
wolffd@0: else % reject the null hypothesis that rho = 0
wolffd@0:   CI = 0;
wolffd@0: end
wolffd@0: p = normcdf(W);
wolffd@0: %p = mynormcdf(W);
wolffd@0: 
wolffd@0: %%%%%%%%%
wolffd@0: 
wolffd@0: function p = normcdf(x,mu,sigma)
wolffd@0: %NORMCDF Normal cumulative distribution function (cdf).
wolffd@0: %   P = NORMCDF(X,MU,SIGMA) computes the normal cdf with mean MU and
wolffd@0: %   standard deviation SIGMA at the values in X.
wolffd@0: %
wolffd@0: %   The size of P is the common size of X, MU and SIGMA. A scalar input  
wolffd@0: %   functions as a constant matrix of the same size as the other inputs.    
wolffd@0: %
wolffd@0: %   Default values for MU and SIGMA are 0 and 1 respectively.
wolffd@0: 
wolffd@0: %   References:
wolffd@0: %      [1]  M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
wolffd@0: %      Functions", Government Printing Office, 1964, 26.2.
wolffd@0: 
wolffd@0: %   Copyright (c) 1993-98 by The MathWorks, Inc.
wolffd@0: %   $Revision: 1.1.1.1 $  $Date: 2005/04/26 02:29:18 $
wolffd@0: 
wolffd@0: if nargin < 3, 
wolffd@0:     sigma = 1;
wolffd@0: end
wolffd@0: 
wolffd@0: if nargin < 2;
wolffd@0:     mu = 0;
wolffd@0: end
wolffd@0: 
wolffd@0: [errorcode x mu sigma] = distchck(3,x,mu,sigma);
wolffd@0: 
wolffd@0: if errorcode > 0
wolffd@0:     error('Requires non-scalar arguments to match in size.');
wolffd@0: end
wolffd@0: 
wolffd@0: %   Initialize P to zero.
wolffd@0: p = zeros(size(x));
wolffd@0: 
wolffd@0: % Return NaN if SIGMA is not positive.
wolffd@0: k1 = find(sigma <= 0);
wolffd@0: if any(k1)
wolffd@0:     tmp   = NaN;
wolffd@0:     p(k1) = tmp(ones(size(k1))); 
wolffd@0: end
wolffd@0: 
wolffd@0: % Express normal CDF in terms of the error function.
wolffd@0: k = find(sigma > 0);
wolffd@0: if any(k)
wolffd@0:     p(k) = 0.5 * erfc( - (x(k) - mu(k)) ./ (sigma(k) * sqrt(2)));
wolffd@0: end
wolffd@0: 
wolffd@0: % Make sure that round-off errors never make P greater than 1.
wolffd@0: k2 = find(p > 1);
wolffd@0: if any(k2)
wolffd@0:     p(k2) = ones(size(k2));
wolffd@0: end
wolffd@0: 
wolffd@0: %%%%%%%%
wolffd@0: 
wolffd@0: function x = norminv(p,mu,sigma);
wolffd@0: %NORMINV Inverse of the normal cumulative distribution function (cdf).
wolffd@0: %   X = NORMINV(P,MU,SIGMA) finds the inverse of the normal cdf with
wolffd@0: %   mean, MU, and standard deviation, SIGMA.
wolffd@0: %
wolffd@0: %   The size of X is the common size of the input arguments. A scalar input  
wolffd@0: %   functions as a constant matrix of the same size as the other inputs.    
wolffd@0: %
wolffd@0: %   Default values for MU and SIGMA are 0 and 1 respectively.
wolffd@0: 
wolffd@0: %   References:
wolffd@0: %      [1]  M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
wolffd@0: %      Functions", Government Printing Office, 1964, 7.1.1 and 26.2.2
wolffd@0: 
wolffd@0: %   Copyright (c) 1993-98 by The MathWorks, Inc.
wolffd@0: %   $Revision: 1.1.1.1 $  $Date: 2005/04/26 02:29:18 $
wolffd@0: 
wolffd@0: if nargin < 3, 
wolffd@0:     sigma = 1;
wolffd@0: end
wolffd@0: 
wolffd@0: if nargin < 2;
wolffd@0:     mu = 0;
wolffd@0: end
wolffd@0: 
wolffd@0: [errorcode p mu sigma] = distchck(3,p,mu,sigma);
wolffd@0: 
wolffd@0: if errorcode > 0
wolffd@0:     error('Requires non-scalar arguments to match in size.');
wolffd@0: end
wolffd@0: 
wolffd@0: % Allocate space for x.
wolffd@0: x = zeros(size(p));
wolffd@0: 
wolffd@0: % Return NaN if the arguments are outside their respective limits.
wolffd@0: k = find(sigma <= 0 | p < 0 | p > 1);
wolffd@0: if any(k)
wolffd@0:     tmp  = NaN;
wolffd@0:     x(k) = tmp(ones(size(k))); 
wolffd@0: end
wolffd@0: 
wolffd@0: % Put in the correct values when P is either 0 or 1.
wolffd@0: k = find(p == 0);
wolffd@0: if any(k)
wolffd@0:     tmp  = Inf;
wolffd@0:     x(k) = -tmp(ones(size(k)));
wolffd@0: end
wolffd@0: 
wolffd@0: k = find(p == 1);
wolffd@0: if any(k)
wolffd@0:     tmp  = Inf;
wolffd@0:     x(k) = tmp(ones(size(k))); 
wolffd@0: end
wolffd@0: 
wolffd@0: % Compute the inverse function for the intermediate values.
wolffd@0: k = find(p > 0  &  p < 1 & sigma > 0);
wolffd@0: if any(k),
wolffd@0:     x(k) = sqrt(2) * sigma(k) .* erfinv(2 * p(k) - 1) + mu(k);
wolffd@0: end