%% (a) bootstrap:
nperm = 1000;
rb = zeros(nperm,1);
for i=1:nperm
    % indices for resampling the data:
    inx = randi(length(x), length(x), 1);
    % resampled data pairs:
    xb=x(inx);
    yb=y(inx);
    rb(i) = corr(xb, yb);
end

%% (b) pdf of the correlation coefficients:
[hb,bb] = hist(rb, 20);
hb = hb/sum(hb)/(bb(2)-bb(1));  % normalization

%% (c) significance:
rbq = quantile(rb, 0.05);
fprintf('correlation coefficient at 5%% significance = %.2f\n', rbq);
if rbq > 0.0
    fprintf('--> correlation r=%.2f is significant\n', rd);
else
    fprintf('--> r=%.2f is not a significant correlation\n', rd);
end

%% plot:
hold on;
bar(b, h, 'facecolor', [0.5 0.5 0.5]);
bar(bb, hb, 'facecolor', 'b');
bar(bb(bb<=rbq), hb(bb<=rbq), 'facecolor', 'r');
plot([rd rd], [0 4], 'r', 'linewidth', 2);
xlim([-0.25 0.75])
xlabel('Correlation coefficient');
ylabel('Probability density');
hold off;

savefigpdf(gcf, 'correlationbootstrap.pdf', 12, 6);