%% Parametric and nonparametric (randomization) test of significance of
% correlation coefficient for simulated data.
m = 5; % Number of variables
n = 100; % Sample size
nullrho = 0; % If 1 there are no cross-correlations (the correlation matrix,
% is unitary), otherwise there are correlations given by the off diagonal
% components of a correlation matrix, randomly drawn.
exptrans = 1; % the power transform of the generated samples.
nsur = 1000; % number of randomized samples.
pridir = 'C:\MyFiles\Teach\Networks\Figures\';
toplot = 0;
toprint = 0;

%% Generate the original sample
if nullrho
    dattxt = sprintf('GaussRhoZeroExp%dn%d',exptrans,m);
else
    dattxt = sprintf('GaussRhoNonZeroExp%dn%d',exptrans,m);
end
if nullrho
    rhoM = eye(m,m);
else
    rhoM = gallery('randcorr',5);
end
muV = zeros(m,1);
xM = mvnrnd(muV,rhoM,n); % Generate multivariate normal sample of size n
xM = xM'.^exptrans; % Power transform of the sample.

%% Declare variables to store results in
count = 0;
rM = NaN*ones(m,m);
statM = NaN*ones(m*(m-1)/2,4); % col 1,2->t-stat,p-value, col 2,3->z-stat,p-value
prnkM = NaN*ones(m*(m-1)/2,2); % col 1-> order, col 2-> p-value
rhoV = NaN*ones(m*(m-1)/2,1); % The true correlation coefficients
rV = NaN*ones(m*(m-1)/2,1); % The estimated correlation coefficients

%% Run for each pair of variables
for i=1:m-1
    for j=i+1:m
        count = count+1;
        %% Compute original r 
        fprintf('pair (%d,%d) ... \n',i,j);
        tmpM = corrcoef(xM(i,:),xM(j,:));
        rM(i,j)=tmpM(1,2);
        rV(count) = rM(i,j);
        rhoV(count) = rhoM(i,j);
        %% Make a z and t parametric test of significance for rho 
        statM(count,1) = sqrt(n-2)*rM(i,j)/sqrt(1-rM(i,j).^2);
        statM(count,2) = 2*(1-tcdf(abs(statM(count,1)),n-2));
        statM(count,3) = sqrt(n-3)*0.5*log((1+rM(i,j))/(1-rM(i,j)));
        statM(count,4) = 2*(1-normcdf(abs(statM(count,3))));
        %% Compute randomized samples and make nonparametric test
        rsurV = NaN*ones(nsur+1,1);
        rsurV(1)=rM(i,j);
        for isur=1:nsur
            yV = xM(j,randperm(n)); % Random suffling of the sample of j variable
            tmpM = corrcoef(xM(i,:),yV);
            rsurV(isur+1)=tmpM(1,2);
        end
        [orV, irV]=sort(rsurV); % Sort the M+1 r-values
        prnkM(count,1)=find(irV==1); % and find the rank of the original
        % Instead of using the simple formula rank/nsur to define the
        % probability being at the tail of the distribution, a more
        % advanced (and presumably more correct) formula is used.
        if prnkM(count,1)>nsur/2
            prnkM(count,2) = 2*(1-(prnkM(count,1)-0.326)/(nsur+1+0.348));
        else
            prnkM(count,2) = 2*(prnkM(count,1)-0.326)/(nsur+1+0.348);
        end
    end
    if toplot
        %% Make a plot of the M+1 r-values for the (i,j) variable pair.
        figure(count)
        clf
        hist(rsurV(2:end));
        hold on
        ax = axis;
        plot(rsurV(1)*[1 1],[ax(3) ax(4)],'r')
        text(ax(1)+0.02*(ax(2)-ax(1)),ax(4)-0.1*(ax(4)-ax(3)),...
            sprintf('r_{%d,%d}',i,j))
        if toprint
            eval(['print -djpeg90 ',pridir,dattxt,'rrand',int2str(count),'.jpg'])
        end
    end
end
%% Display the results of the parametric and nonparametric test
fprintf('rho \t r \n');
disp([rhoV rV])
fprintf('t-stat \t z-stat \t rank \n');
disp([statM(:,[1 3]) prnkM(:,1)])
fprintf('t-stat p-val \t z-stat p-val \t rank p-val \n');
disp([statM(:,[2 4]) prnkM(:,2)])





