-
Notifications
You must be signed in to change notification settings - Fork 1
/
pub_atten_amplitude.py
executable file
·124 lines (96 loc) · 4.12 KB
/
pub_atten_amplitude.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
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
#!/usr/bin/python
"""
Plot best-fit
"""
import numpy as np
import pylab as P
import scipy.interpolate
import scipy.integrate
import scipy.optimize
# Peak wavelengths (ugriz) in Angstrom
l = [3543., 4770., 6231., 7625., 9134.]
BANDS = ['u', 'g', 'r', 'i', 'z']
colours = ['#DA49DA', '#00a5ff', '#ff8200', '#F34334', '#9D2020']
# Data from scan_attenuation.py
# Best-fit, median, 16% CL, 84% CL, -1sig, +1sig
#dat_u = [5.1596e+01, 5.264e+01, 4.654e+01, 5.976e+01, 6.108e+00, 7.119e+00]
## u-band best-fit amp: 5.1596e+01 (logL = -75.5)
#dat_g = [5.8313e-01, 5.656e-01, 4.699e-01, 6.605e-01, 9.572e-02, 9.489e-02]
## g-band best-fit amp: 5.8313e-01 (logL = -23.8)
#dat_r = [9.0460e-03, 9.229e-03, 7.357e-03, 1.121e-02, 1.872e-03, 1.981e-03]
## r-band best-fit amp: 9.0460e-03 (logL = -80.3)
#dat_i = [3.2924e-06, 3.930e-06, 2.416e-06, 6.080e-06, 1.514e-06, 2.150e-06]
## i-band best-fit amp: 3.2924e-06 (logL = -20.8)
#dat_z = [8.7386e-07, 1.421e-06, 7.015e-07, 2.436e-06, 7.199e-07, 1.015e-06]
## z-band best-fit amp: 8.7386e-07 (logL = -21.4)
#######################
# Load likelihood data from saved files
logl_interp = []
median = []; errp = []; errm = []
for b in BANDS:
# Load cache likelihood data generated by scan_attenuation.py
pvals, logl = np.load("atten_like_%s.npy" % b).T
# Interpolate the likelihood and store interpolation fn. in list
_logl_interp = scipy.interpolate.interp1d(np.log(pvals), logl,
kind='quadratic', bounds_error=False,
fill_value=-500.)
pp = np.logspace(np.log10(np.min(pvals)), np.log10(np.max(pvals)), 1000)
_logl = _logl_interp(np.log(pp))
logl_interp.append(_logl_interp)
# Calculate cumulative dist. fn.
cdf = scipy.integrate.cumtrapz(np.exp(_logl - np.max(_logl)), pp, initial=0.)
cdf /= cdf[-1]
p_cdf = scipy.interpolate.interp1d(cdf, pp, kind='linear')
median.append(p_cdf(0.5))
errp.append(p_cdf(0.84)- p_cdf(0.5))
errm.append(p_cdf(0.5)- p_cdf(0.16))
def extinct_amp(l, params):
"""
Simple model for the extinction amplitude as a function of wavelength.
"""
amp, scale, lambda0 = params
return amp * np.exp(-scale * (l - lambda0))
# Try to fit function of frequency, using stored likelihoods
def like(params):
amp, scale, lambda0 = params
logl = []
for i in range(len(BANDS)):
# Get interpolated log-like for this extinction amplitude in this band
_extinct_amp = extinct_amp(l[i], params)
logl.append( logl_interp[i](np.log(_extinct_amp)) )
return logl
params0 = [0.5, 0.003, 5000.]
pnew = scipy.optimize.leastsq(like, params0)[0]
for pname, bfval in zip(['tau0', 'kappa', 'lambda0'], pnew):
print "%10s: %f" % (pname, bfval)
print "logL(best-fit) = %3.1f" % np.sum(like(pnew))
## Calculate median and 1/2-sigma bounds
#print "%3.3e %3.3e %3.3e %3.3e %3.3e" % (p_cdf(0.5), p_cdf(0.16), p_cdf(0.84),
# p_cdf(0.5)- p_cdf(0.16), p_cdf(0.84)- p_cdf(0.5))
#######################
#data = np.array([dat_u, dat_g, dat_r, dat_i, dat_z]).T
#print data.shape
P.subplot(111)
#P.plot(l, data[0], 'kx', ms=8., mew=1.5)
P.errorbar(l, median, yerr=(errm, errp), color='#0B73CE', ms=8.,
ls='none', marker='.', capsize=4., elinewidth=1.5, mew=1.5)
# Plot best-fit powerlaw
ll = np.linspace(2700., 9900., 100)
#P.plot(ll, 3e6*np.exp(-16.*ll/5000.), 'k--', lw=2., alpha=0.5)
#P.plot(ll, 0.5*np.exp(-0.003*(ll - 5000.)), 'k--', lw=2., alpha=0.5)
P.plot(ll, extinct_amp(ll, pnew), color='#E72327', lw=2., dashes=[4,3])
# Place band names close to markers
for i in range(len(BANDS)):
P.annotate(BANDS[i], xy=(l[i], 5e-8), fontsize=22., color=colours[i],
horizontalalignment='center')
P.yscale('log')
P.xlabel(r"$\lambda\, [\AA]$", fontsize=18)
P.ylabel(r"$f\,(\lambda)$", fontsize=18)
P.gca().tick_params(axis='both', which='major', labelsize=20, size=8.,
width=1.5, pad=8.)
P.gca().tick_params(axis='both', which='minor', labelsize=20, size=5.,
width=1.5, pad=8.)
P.xlim((2700., 9700.))
P.tight_layout()
P.savefig('../draft/atten_freq.pdf')
P.show()