In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import pymc
import astropy
import astropy.table as table
import pymc.graph
import SPaCT
from IPython.display import Image
from triangle import corner
plt.rc('text', usetex = False)
objname = 'NGC2558'
fiberfits = table.Table.read('{}/fiberfits.dat'.format(objname), format = 'ascii')
fiberfits = fiberfits[np.isfinite(fiberfits['V'])]
ra = fiberfits['ra']
dec = fiberfits['dec']
V = fiberfits['V']
dV = fiberfits['dV']
sig = fiberfits['sigma']
dsig = fiberfits['dsigma']
sig_inst = 250.
print fiberfits[[0]]
print V[0]
In [49]:
def make_tilted_disk_model(phi_0, i, V_sys, v_rot, h_rot, sig):
ar_len = 40.
disk_d = 60.
RA = np.linspace(-40., 40., 11)
DEC = np.linspace(-40., 40., 11)
RR, DD = np.meshgrid(RA, DEC)
reload(SPaCT)
vmap_mock = SPaCT.tilted_disk_model_tanh_(ra = RR, dec = DD, phi_0 = phi_0, i = i,
V_sys = V_sys, v_rot = v_rot, h_rot = h_rot)
vmap_err = np.random.normal(0., sig, RR.shape)
vmap_mock += vmap_err
#print vmap_err
plt.figure(figsize = (6, 6))
ax = plt.subplot(111)
vm = ax.scatter(RR, DD, c = vmap_mock, cmap = 'RdBu_r', vmin = V_sys - 400, vmax = V_sys + 400,
marker = 'o', edgecolor = 'k', s = 75)
#ax.arrow(0., 0., ar_len*np.cos(phi_0), ar_len*np.sin(phi_0), zorder = 3,
# head_length = 5., head_width = 5., fc = 'k', ec = 'k') # major axis
if i == 0:
minaxisarrcolor = 'k'
elif i < 0:
minaxisarrcolor = 'r'
else:
minaxisarrcolor = 'b'
#ax.arrow(0., 0., ar_len*np.cos(phi_0 + np.pi/2.)*np.cos(i), ar_len*np.sin(phi_0 + np.pi/2.)*np.cos(i), zorder = 3,
# head_length = 5. + 2.5*np.sin(i), head_width = 5. + 2.5*np.sin(i),
# fc = minaxisarrcolor, ec = minaxisarrcolor) # minor axis
from matplotlib.patches import Ellipse
e = Ellipse(xy = (0,0), width = disk_d, height = disk_d*np.cos(i),
angle = np.degrees(phi_0), zorder = 2, fc = 'None', ec = 'k', linewidth = 5)
#ax.add_artist(e)
e = Ellipse(xy = (0,0), width = 0.5*disk_d, height = 0.5*disk_d*np.cos(i),
angle = np.degrees(phi_0), zorder = 2, fc = 'None', ec = 'grey', linewidth = 3)
#ax.add_artist(e)
ax.text(-50., -50., r'$v_{rot}=$' + format(v_rot) + ' km/s', size = 14)
ax.text(5., -50., r'$h_{rot}=$' + format(h_rot) + ' arcsec', size = 14)
ax.set_aspect('equal')
ax.set_xlim([-55., 55.])
ax.set_ylim([-55., 55.])
ax.set_title(r'Tilted disk: $i=${0:.2f} $\pi$, $\phi_0=${1:.2f} $\pi$'.format(i/np.pi, phi_0/np.pi))
ax.set_xlabel('RA offset [arcsec]')
ax.set_ylabel('Dec offset [arcsec]')
cb = plt.colorbar(vm, shrink = 0.8)
cb.set_label('$V_{tot}$ [km/s]', size = 16)
plt.show()
return vmap_mock
phi_0_actual = np.pi/4.
incl_actual = np.pi/3.
V_sys_actual = 0.
v_rot_actual = 400.
h_rot_actual = 20.
vsig_actual = 0.2*sig_inst
vsig_actual = 25.
vmap_mock = make_tilted_disk_model(phi_0 = phi_0_actual, i = incl_actual, V_sys = V_sys_actual,
v_rot = v_rot_actual, h_rot = h_rot_actual, sig = vsig_actual)
vmap_mock = make_tilted_disk_model(phi_0 = phi_0_actual, i = incl_actual, V_sys = V_sys_actual,
v_rot = v_rot_actual, h_rot = 2.*h_rot_actual, sig = vsig_actual)
vmap_mock = make_tilted_disk_model(phi_0 = phi_0_actual, i = 0, V_sys = V_sys_actual,
v_rot = v_rot_actual, h_rot = h_rot_actual, sig = vsig_actual)
vmap_mock = make_tilted_disk_model(phi_0 = phi_0_actual, i = np.pi/2, V_sys = V_sys_actual,
v_rot = v_rot_actual, h_rot = h_rot_actual, sig = vsig_actual)
vmap_mock = make_tilted_disk_model(phi_0 = 0., i = np.pi/2, V_sys = V_sys_actual,
v_rot = v_rot_actual, h_rot = h_rot_actual, sig = vsig_actual)
vmap_mock = make_tilted_disk_model(phi_0 = np.pi/2, i = np.pi/2, V_sys = V_sys_actual,
v_rot = v_rot_actual, h_rot = h_rot_actual, sig = vsig_actual)
In [13]:
#hyperpriors
#v_rot_mu = pymc.TruncatedNormal('v_rot_mu', mu = 300., tau = 1./500**2., a = 0., b = 1000.)
#v_rot_sig = pymc.TruncatedNormal('v_rot_sig', mu = 300., tau = 1./500**2., a = 0., b = 1000.)
#priors
'''
the periodicity of priors is an issue, so uniform priors from 0 to 2pi can be transformed by chained
cos and arccos operations on the chains themselves.
'''
phi_0 = pymc.Uniform('phi_0', 0., 2.*np.pi)
incl = pymc.Uniform('incl', 0., 2.*np.pi)
vmap_mock_ctr = vmap_mock[vmap_mock.shape[0]/2, vmap_mock.shape[1]/2]
V_sys = pymc.Normal('V_sys', mu = vmap_mock_ctr,
tau = 1./vsig_actual**2., value = vmap_mock_ctr)
v_obs_err = pymc.HalfNormal('v_obs_err', tau = 1./vsig_actual**2.)
#v_rot = pymc.TruncatedNormal('v_rot', mu = v_rot_mu, tau = 1./v_rot_sig**2., a = 0.)
v_rot = pymc.TruncatedNormal('v_rot', mu = 300., tau = 1./500.**2., a = 0., b = 700.)
h_rot = pymc.TruncatedNormal('h_rot', mu = 10., tau = 1./50**2., a = 0.)
#model
@pymc.deterministic
def tilted_disk_model_tanh(ra = ra, dec = dec, phi_0 = phi_0, incl = incl, V_sys = V_sys, v_rot = v_rot, h_rot = h_rot):
reload(SPaCT)
return SPaCT.tilted_disk_model_tanh_(ra, dec, phi_0, incl, V_sys, v_rot, h_rot)
#likelihood: Normal distribution, of stdev dV
V_likelihood = pymc.Normal('V', tilted_disk_model_tanh, 1./v_obs_err**2., value = V, observed = True)
#posterior predictive
V_pred = pymc.Normal('V', tilted_disk_model_tanh, 1./v_obs_err**2.)
model_tilted_disk_tanh = pymc.Model({'incl': incl, 'phi_0': phi_0, 'V_sys': V_sys, 'v_rot': v_rot,
'V_likelihood': V_likelihood, 'V_pred': V_pred, 'v_obs_err': v_obs_err})
graph = pymc.graph.graph(model_tilted_disk_tanh)
graph.write_png('tilted_disk_tanh_model.png')
mcmc = pymc.MCMC(model_tilted_disk_tanh)
Image('tilted_disk_tanh_model.png')
Out[13]:
In [14]:
nchains = 1
for _ in range(nchains):
mcmc.sample(iter = 10000, burn = 1000)
In [15]:
pymc.Matplot.plot(mcmc)
In [30]:
corner_incl = ['incl', 'phi_0', 'v_obs_err', 'h_rot', 'v_rot', 'V_sys']
angular = [True, True, False, False, False, False]
truths = [incl_actual, phi_0_actual, vsig_actual, h_rot_actual, v_rot_actual, V_sys_actual]
samples = np.array([mcmc.trace(label, chain=None)[:] for label in corner_incl]).T
for (i, angular) in enumerate(samples.T):
#print i
if angular[i] == True:
samples[i] = np.arccos(np.cos(samples.T[i]))
plt.figure(figsize = (10, 10))
corner(samples, labels = corner_incl, quantiles=[0.16, 0.5, 0.84], truths = truths)
plt.show()
for i, row in enumerate(samples.T):
print '{0}: {1}'.format(corner_incl[i], truths[i])
print np.percentile(row, [16, 50, 84])
In [ ]: