/
cos_pymcc.py
83 lines (73 loc) · 2.73 KB
/
cos_pymcc.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
#!/Users/drupke/opt/anaconda3/envs/astroconda/bin python
from astropy.io import ascii
from astropy.table import Table
import numpy as np
from pymccorrelation import pymccorrelation
Nperturb = 1000
datadir = '/Users/drupke/Box Sync/cosquest/plots/correlations/'
data = ascii.read(datadir+'weq_vs_nhxray_dat.txt')
res = pymccorrelation(data['xdat'],data['ydat'],data['xerr'],data['yerr'],
data['xl'],data['yl'],coeff='kendallt',Nperturb=Nperturb,
return_dist=True)
# Compute pval from cc distribution
scc = np.sort(res[2])
izero = np.searchsorted(scc,0.)
# case of positive cc, all > 0
# case of negative cc, all < 0
if izero == 0 or izero == scc.size-1:
pvaldist = 1./float(Nperturb)
pvalul = True
else:
pvalul = False
pvaldist = float(izero)/float(scc.size)
if res[0][1] < 0.:
pvaldist = 1. - pvaldist
outtab = Table()
outtab['cc:16/50/84'] = res[0]
outtab['pval:16/50/84'] = res[1]
outtab['pval:calc/ul'] = [pvaldist,pvalul,'N/A']
outtab.write(datadir+'weq_vs_nhxray_stat.txt', format='ascii', delimiter='\t', overwrite=True)
data = ascii.read(datadir+'vwtavg_vs_nhxray_dat.txt')
res = pymccorrelation(data['xdat'],data['ydat'],data['xerr'],data['yerr'],
data['xl'],data['yl'],coeff='kendallt',Nperturb=Nperturb,
return_dist=True)
# Compute pval from cc distribution
scc = np.sort(res[2])
izero = np.searchsorted(scc,0.)
# case of positive cc, all > 0
# case of negative cc, all < 0
if izero == 0 or izero == scc.size-1:
pvaldist = 1./float(Nperturb)
pvalul = True
else:
pvalul = False
pvaldist = float(izero)/float(scc.size)
if res[0][1] < 0.:
pvaldist = 1. - pvaldist
outtab = Table()
outtab['cc:16/50/84'] = res[0]
outtab['pval:16/50/84'] = res[1]
outtab['pval:calc/ul'] = [pvaldist,pvalul,'N/A']
outtab.write(datadir+'vwtavg_vs_nhxray_stat.txt', format='ascii', delimiter='\t', overwrite=True)
data = ascii.read(datadir+'vwtrms_vs_nhxray_dat.txt')
res = pymccorrelation(data['xdat'],data['ydat'],data['xerr'],data['yerr'],
data['xl'],data['yl'],coeff='kendallt',Nperturb=Nperturb,
return_dist=True)
# Compute pval from cc distribution
scc = np.sort(res[2])
izero = np.searchsorted(scc,0.)
# case of positive cc, all > 0
# case of negative cc, all < 0
if izero == 0 or izero == scc.size-1:
pvaldist = 1./float(Nperturb)
pvalul = True
else:
pvalul = False
pvaldist = float(izero)/float(scc.size)
if res[0][1] < 0.:
pvaldist = 1. - pvaldist
outtab = Table()
outtab['cc:16/50/84'] = res[0]
outtab['pval:16/50/84'] = res[1]
outtab['pval:calc/ul'] = [pvaldist,pvalul,'N/A']
outtab.write(datadir+'vwtrms_vs_nhxray_stat.txt', format='ascii', delimiter='\t', overwrite=True)