91 lines
3.1 KiB
Python
91 lines
3.1 KiB
Python
import numpy as np
|
|
import scipy.stats as st
|
|
import matplotlib.pyplot as plt
|
|
import matplotlib.gridspec as gridspec
|
|
import matplotlib.ticker as ticker
|
|
from plotstyle import *
|
|
|
|
rng = np.random.RandomState(37281)
|
|
|
|
# generate correlated data:
|
|
n = 200
|
|
a = 0.2
|
|
x = rng.randn(n);
|
|
y = rng.randn(n) + a*x;
|
|
#x = rng.exponential(1.0, n);
|
|
#y = rng.exponential(2.0, n) + a*x;
|
|
|
|
rd = np.corrcoef(x, y)[0, 1]
|
|
|
|
# permutation:
|
|
nperm = 1000
|
|
rs = []
|
|
for i in range(nperm) :
|
|
xr = rng.permutation(x)
|
|
yr = rng.permutation(y)
|
|
rs.append(np.corrcoef(xr, yr)[0, 1])
|
|
rr = np.corrcoef(xr, yr)[0, 1]
|
|
|
|
# pdf of the correlation coefficients:
|
|
h, b = np.histogram(rs, 20, density=True)
|
|
|
|
# significance:
|
|
rq = np.percentile(rs, 95.0)
|
|
print('Measured correlation coefficient = %.2f, correlation coefficient at 95%% percentile of permutation = %.2f' % (rd, rq))
|
|
ra = 1.0-0.01*st.percentileofscore(rs, rd)
|
|
print('Measured correlation coefficient %.2f is at %.4f percentile of permutation' % (rd, ra))
|
|
|
|
rp, ra = st.pearsonr(x, y)
|
|
print('Measured correlation coefficient %.2f is at %.4f percentile of test' % (rp, ra))
|
|
|
|
fig = plt.figure(figsize=cm_size(figure_width, 1.8*figure_height))
|
|
gs = gridspec.GridSpec(nrows=2, ncols=2, wspace=0.35, hspace=0.8,
|
|
**adjust_fs(fig, left=5.0, right=1.5, top=1.0, bottom=2.7))
|
|
|
|
ax = fig.add_subplot(gs[0,0])
|
|
ax.text(0.0, 4.0, 'r=%.2f' % rd, ha='center')
|
|
ax.plot(x, y, **psAm)
|
|
ax.set_xlim(-4.0, 4.0)
|
|
ax.set_ylim(-4.0, 4.0)
|
|
ax.xaxis.set_major_locator(ticker.MultipleLocator(2.0))
|
|
ax.yaxis.set_major_locator(ticker.MultipleLocator(2.0))
|
|
ax.set_xlabel('x')
|
|
ax.set_ylabel('y')
|
|
|
|
ax = fig.add_subplot(gs[0,1])
|
|
ax.text(0.0, 4.0, 'r=%.2f' % rr, ha='center')
|
|
ax.plot(xr, yr, **psAm)
|
|
ax.set_xlim(-4.0, 4.0)
|
|
ax.set_ylim(-4.0, 4.0)
|
|
ax.xaxis.set_major_locator(ticker.MultipleLocator(2.0))
|
|
ax.yaxis.set_major_locator(ticker.MultipleLocator(2.0))
|
|
ax.set_xlabel('Shuffled x')
|
|
ax.set_ylabel('Shuffled y')
|
|
|
|
ax = fig.add_subplot(gs[1,:])
|
|
ax.annotate('Measured\ncorrelation\nis significant!',
|
|
xy=(rd, 1.1), xycoords='data',
|
|
xytext=(rd-0.01, 3.0), textcoords='data', ha='left',
|
|
arrowprops=dict(arrowstyle="->", relpos=(0.3,0.0),
|
|
connectionstyle="angle3,angleA=10,angleB=80") )
|
|
ax.annotate('95% percentile',
|
|
xy=(0.14, 0.9), xycoords='data',
|
|
xytext=(0.16, 6.2), textcoords='data', ha='left',
|
|
arrowprops=dict(arrowstyle="->", relpos=(0.1,0.0),
|
|
connectionstyle="angle3,angleA=30,angleB=80") )
|
|
ax.annotate('Distribution of\nuncorrelated\nsamples',
|
|
xy=(-0.08, 3.6), xycoords='data',
|
|
xytext=(-0.22, 5.0), textcoords='data', ha='left',
|
|
arrowprops=dict(arrowstyle="->", relpos=(0.5,0.0),
|
|
connectionstyle="angle3,angleA=150,angleB=110") )
|
|
ax.bar(b[:-1], h, width=b[1]-b[0], **fsC)
|
|
ax.bar(b[:-1][b[:-1]>=rq], h[b[:-1]>=rq], width=b[1]-b[0], **fsB)
|
|
ax.plot( [rd, rd], [0, 1], **lsA)
|
|
ax.set_xlim(-0.25, 0.35)
|
|
ax.set_ylim(0.0, 6.0)
|
|
ax.yaxis.set_major_locator(ticker.MultipleLocator(2.0))
|
|
ax.set_xlabel('Correlation coefficient')
|
|
ax.set_ylabel('PDF of H0')
|
|
|
|
plt.savefig('permutecorrelation.pdf')
|