function [CI, r, p] = cond_indep_fisher_z(X, Y, S, C, N, alpha)
% COND_INDEP_FISHER_Z Test if X indep Y given Z using Fisher's Z test
% CI = cond_indep_fisher_z(X, Y, S, C, N, alpha)
%
% C is the covariance (or correlation) matrix
% N is the sample size
% alpha is the significance level (default: 0.05)
%
% See p133 of T. Anderson, "An Intro. to Multivariate Statistical Analysis", 1984

if nargin < 6, alpha = 0.05; end

r = partial_corr_coef(C, X, Y, S);
z = 0.5*log( (1+r)/(1-r) );
z0 = 0;
W = sqrt(N - length(S) - 3)*(z-z0); % W ~ N(0,1)
cutoff = norminv(1 - 0.5*alpha); % P(|W| <= cutoff) = 0.95
%cutoff = mynorminv(1 - 0.5*alpha); % P(|W| <= cutoff) = 0.95
if abs(W) < cutoff
  CI = 1;
else % reject the null hypothesis that rho = 0
  CI = 0;
end
p = normcdf(W);
%p = mynormcdf(W);

%%%%%%%%%

function p = normcdf(x,mu,sigma)
%NORMCDF Normal cumulative distribution function (cdf).
%   P = NORMCDF(X,MU,SIGMA) computes the normal cdf with mean MU and
%   standard deviation SIGMA at the values in X.
%
%   The size of P is the common size of X, MU and SIGMA. A scalar input  
%   functions as a constant matrix of the same size as the other inputs.    
%
%   Default values for MU and SIGMA are 0 and 1 respectively.

%   References:
%      [1]  M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
%      Functions", Government Printing Office, 1964, 26.2.

%   Copyright (c) 1993-98 by The MathWorks, Inc.
%   $Revision: 1.1 $  $Date: 2004/02/10 18:58:56 $

if nargin < 3, 
    sigma = 1;
end

if nargin < 2;
    mu = 0;
end

[errorcode x mu sigma] = distchck(3,x,mu,sigma);

if errorcode > 0
    error('Requires non-scalar arguments to match in size.');
end

%   Initialize P to zero.
p = zeros(size(x));

% Return NaN if SIGMA is not positive.
k1 = find(sigma <= 0);
if any(k1)
    tmp   = NaN;
    p(k1) = tmp(ones(size(k1))); 
end

% Express normal CDF in terms of the error function.
k = find(sigma > 0);
if any(k)
    p(k) = 0.5 * erfc( - (x(k) - mu(k)) ./ (sigma(k) * sqrt(2)));
end

% Make sure that round-off errors never make P greater than 1.
k2 = find(p > 1);
if any(k2)
    p(k2) = ones(size(k2));
end

%%%%%%%%

function x = norminv(p,mu,sigma);
%NORMINV Inverse of the normal cumulative distribution function (cdf).
%   X = NORMINV(P,MU,SIGMA) finds the inverse of the normal cdf with
%   mean, MU, and standard deviation, SIGMA.
%
%   The size of X is the common size of the input arguments. A scalar input  
%   functions as a constant matrix of the same size as the other inputs.    
%
%   Default values for MU and SIGMA are 0 and 1 respectively.

%   References:
%      [1]  M. Abramowitz and I. A. Stegun, "Handbook of Mathematical
%      Functions", Government Printing Office, 1964, 7.1.1 and 26.2.2

%   Copyright (c) 1993-98 by The MathWorks, Inc.
%   $Revision: 1.1 $  $Date: 2004/02/10 18:58:56 $

if nargin < 3, 
    sigma = 1;
end

if nargin < 2;
    mu = 0;
end

[errorcode p mu sigma] = distchck(3,p,mu,sigma);

if errorcode > 0
    error('Requires non-scalar arguments to match in size.');
end

% Allocate space for x.
x = zeros(size(p));

% Return NaN if the arguments are outside their respective limits.
k = find(sigma <= 0 | p < 0 | p > 1);
if any(k)
    tmp  = NaN;
    x(k) = tmp(ones(size(k))); 
end

% Put in the correct values when P is either 0 or 1.
k = find(p == 0);
if any(k)
    tmp  = Inf;
    x(k) = -tmp(ones(size(k)));
end

k = find(p == 1);
if any(k)
    tmp  = Inf;
    x(k) = tmp(ones(size(k))); 
end

% Compute the inverse function for the intermediate values.
k = find(p > 0  &  p < 1 & sigma > 0);
if any(k),
    x(k) = sqrt(2) * sigma(k) .* erfinv(2 * p(k) - 1) + mu(k);
end