function [CorrAvg, pValue] = qapBayesTest(IONetwork,EqNetwDraws,QAPReps)

CorrAvg = 0;
pValue = 0;
for ii = 1:1:QAPReps
    %Select a random Gibbs draw
        EqNetw = squeeze(EqNetwDraws(:,:,randi(size(EqNetwDraws,3))));
    %Actual correl with draw
        ActCorrWithDraw = corr(EqNetw(:), IONetwork(:));
        CorrAvg = CorrAvg + ActCorrWithDraw;
    %Get a random reordering of the IO matrix
        RandOrd = randperm(size(IONetwork,1));
        IOrand = IONetwork(RandOrd,RandOrd);
    %Update pvalue if correlation with actual is lower than with IOrand
        if ActCorrWithDraw < corr(EqNetw(:), IOrand(:))
            pValue = pValue+1;
        end
end

CorrAvg = CorrAvg / QAPReps;
pValue  = pValue  / QAPReps;
end