We approximate the local group distance, radial velocity and proper motion likelihood function by sampling from the posterior distributions for these variables reported in the literature and transforming to kinematic variables in the M31-centric coordinate system.
In [1]:
%matplotlib inline
import localgroup
import triangle
import sklearn
from sklearn import mixture
import numpy as np
import pickle
import matplotlib.patches as mpatches
Inside the Likelihood object is a "triplet" object called T, which contains an array of sample local groups, each with kinematic parameters consistent with the observational data. Let's plot these kinematic parameters in a "triangle" figure, to show all their 1 and 2-D marginal distributions.
In [2]:
L = localgroup.Likelihood(isPair=False)
L.generate(Nsamples=200000)
In [3]:
L.set_PDF(mixture.GMM(n_components=10, covariance_type='full'))
L.approximate()
In [4]:
figure_obs = L.plot_samples(10, color='b', overlay=False)
The above plot shows a Gaussian Mixture model fitted Gaussians. The shaded regions show two standard deviations. The samples data has been preprocessed to zero the mean and scale by standard deviation. Since we are using the Gaussian Mixture Model to model the underlying PDF of the data, more components is always better.
Due to lack of a standard goodness of fit test for GMM's, the best we can do is graphically show that the model reproduces the data well. We proceed by drawing a set of points from the fitted model, where each point is a local group with (MW_D, MW_vr, MW_vt, M33_D, M33_vr, M33_vt). We then plot the 1D and 2D marginalizations of the drawn point set and show that the marginalizations match the marginalizations of the true data.
In [5]:
figure_model = L.model_gof(L.T.Nsamples, color="r", fig=None)
In [6]:
figure_model
Out[6]:
In [7]:
L.model_gof(L.T.Nsamples, color="r", fig=figure_obs)
red_patch = mpatches.Patch(color='red')
blue_patch = mpatches.Patch(color='blue')
figure_obs.legend(handles=[red_patch, blue_patch], labels=["Model Generated", "Observation Generated"], fontsize=16)
Out[7]:
In [8]:
figure_obs.suptitle("Likelihood Model Goodness of Fit", fontsize=16)
Out[8]:
In [9]:
figure_obs
Out[9]:
The above plot shows that the points drawn from the model create a population that is very similar to the true data.
In [10]:
#figure_obs.savefig("/afs/slac.stanford.edu/u/ki/mwillia1/Thesis/LocalGroupHaloProps/doc/thesis/plots/final_plots/Trip_LGMM_GOF.pdf", dpi=1200)
Below we read the preconfigured files containing the Consuelo (soon to be Dark Sky) Local Group analogs into a Triplet object. We plot the marginalizations of the simulation data, which allows us to compare with the LG prior.
In [11]:
path = '/lustre/ki/pfs/mwillia1/LG_project/Consuelo_Boxes/All_Boxes_quad_dat_M31_larger.npy'
#path = '/afs/slac.stanford.edu/u/ki/mwillia1/Thesis/data_files/complete_triplets.txt'
#path_054 = '/lustre/ki/pfs/mwillia1/LG_project/Consuelo_Boxes/All_Boxes_trip_dat_M31_larger_0.949.npy'
npoints = 200000
halo_props = ['MW_Mvir', 'M31_Mvir', 'M33_Mvir']
In [12]:
h=np.load(path)
In [13]:
h.size
Out[13]:
In [12]:
Tr = localgroup.Triplet(isPair=False)
#Tr_054 = localgroup.Triplet(isPair=True)
Tr.read_sim_points(path, npoints, halo_props, h=0.7, a=1.0, npy=True)
#Tr_054.read_sim_points(path_054, npoints, halo_props, h=0.7, a=0.949, npy=True)
In [13]:
Tr.transform_to_M31(sim=True)
#Tr_054.transform_to_M31(sim=True)
In [25]:
Tr.sim_samples = np.vstack((Tr.sim_samples, Tr_054.sim_samples))
Tr.MW.Mvir = np.hstack((Tr.MW.Mvir, Tr_054.MW.Mvir))
Tr.M31.Mvir = np.hstack((Tr.M31.Mvir, Tr_054.M31.Mvir))
#Tr.M33.Mvir = np.hstack((Tr.M33.Mvir, Tr_054.M33.Mvir))
Tr.MW.Rvir = np.hstack((Tr.MW.Rvir, Tr_054.MW.Rvir))
Tr.M31.Rvir = np.hstack((Tr.M31.Rvir, Tr_054.M31.Rvir))
#Tr.M33.Rvir = np.hstack((Tr.M33.Rvir, Tr_054.M33.Rvir))
Tr.MW.Rs = np.hstack((Tr.MW.Rs, Tr_054.MW.Rs))
Tr.M31.Rs = np.hstack((Tr.M31.Rs, Tr_054.M31.Rs))
#Tr.M33.Rs = np.hstack((Tr.M33.Rs, Tr_054.M33.Rs))
Tr.MW.Cvir = np.hstack((Tr.MW.Cvir, Tr_054.MW.Cvir))
Tr.M31.Cvir = np.hstack((Tr.M31.Cvir, Tr_054.M31.Cvir))
#Tr.M33.Cvir = np.hstack((Tr.M33.Cvir, Tr_054.M33.Cvir))
Tr.LG_Mvir = np.hstack((Tr.LG_Mvir, Tr_054.LG_Mvir))
In [14]:
Tr.sim_samples.shape
Out[14]:
In [17]:
#Tr.mass_filter('sim')
In [15]:
Tr.dist_filter((Tr.sim_samples[:,0] < 10) & (Tr.sim_samples[:,3] < 10) & (Tr.sim_samples[:,6] < 10))
#Tr.dist_filter((Tr.sim_samples[:,0] < 10))
Out[15]:
In [29]:
max(Tr.sim_samples[:,3])
In [16]:
Tr.preprocess(L.samples_means, L.samples_stds, mode='sim')
In [17]:
sim_plot = Tr.plot_kinematics('sim', L.samples_means, L.samples_stds, color='c', fig=None)
In [ ]:
In [18]:
Tr.plot_kinematics('sim', L.samples_means, L.samples_stds, color='g', fig=figure_model)
Out[18]:
In [19]:
red_patch = mpatches.Patch(color='red')
green_patch = mpatches.Patch(color='green')
figure_model.legend(handles=[red_patch, green_patch], labels=["Model Generated", "Consuelo Triplets"], fontsize=16)
figure_model.suptitle("Model Likelihood Overlayed with Consuelo Prior", fontsize=16)
figure_model
Out[19]:
In [ ]:
In [29]:
ax=np.reshape(np.array(figure_model.axes), (9,9))
In [31]:
AX=ax[8,0]
In [33]:
AX.set_xlim((0,10))
Out[33]:
In [34]:
figure_model
Out[34]:
In [19]:
#sim_plot.savefig('/afs/slac.stanford.edu/u/ki/mwillia1/Thesis/LocalGroupHaloProps/doc/thesis/plots/asurps/likelihood_prior.png', dpi=1200)
In [20]:
Tr.compute_model_weights(L, "sim")
In [21]:
N95 = Tr.calculate_N95()
N95
Out[21]:
In [11]:
Tr.weights[Tr.weights[:]>.024797]
Out[11]:
In [15]:
MW_C200 = Tr.calculate_cvir_to_c200(Tr.MW.Cvir)
M31_C200 = Tr.calculate_cvir_to_c200(Tr.M31.Cvir)
M33_C200 = Tr.calculate_cvir_to_c200(Tr.M33.Cvir)
MW_M200 = Tr.calculate_Mvir_to_M200(Tr.MW.Mvir, Tr.MW.Cvir, MW_C200)
M31_M200 = Tr.calculate_Mvir_to_M200(Tr.M31.Mvir, Tr.M31.Cvir, M31_C200)
M33_M200 = Tr.calculate_Mvir_to_M200(Tr.M33.Mvir, Tr.M33.Cvir, M33_C200)
total_200 = M31_M200+MW_M200
In [20]:
#m33 = np.ones(len(Tr.MW.Mvir))
#m33[0]=1.5
#m33[1]=.5
total = Tr.MW.Mvir + Tr.M31.Mvir
all_mvir = np.log10(np.transpose(np.vstack((Tr.MW.Mvir, Tr.M31.Mvir, Tr.M33.Mvir, total))))
#all_mvir = np.log10(np.transpose(np.vstack((MW_M200, M31_M200, M33_M200, total_200))))
#labs = ["$M200_{MW}$", "$M200_{M31}$", "$M200_{M33}$", "$M_{MW+M31}$"]
labs = ["$Mvir_{MW}$", "$Mvir_{M31}$", "$Mvir_{M33}$", "$M_{MW+M31}$"]
In [56]:
all_mvir = np.transpose(np.vstack((Tr.MW.Cvir, Tr.M31.Cvir, Tr.M33.Cvir)))
labs = ["$Cvir_{MW}$", "$Cvir_{M31}$", "$Cvir_{M33}$"]
In [19]:
all_mvir.shape
Out[19]:
In [20]:
figure2 = triangle.corner(all_mvir, labels=labs, quantiles=[0.16,0.5,0.84], fig=None, weights=None,\
plot_contours=True, show_titles=False, title_args={"fontsize": 16}, label_args={"fontsize": 16},\
plot_datapoints=False, bins=20, color='c', alpha=0.1)
#figure.suptitle("Weighted Mass Posterior PDF", fontsize=16, horizontalalignment='left')
In [25]:
figure2
Out[25]:
In [28]:
figure2.savefig('/afs/slac.stanford.edu/u/ki/mwillia1/Thesis/LocalGroupHaloProps/doc/thesis/plots/final_plots/W_vs_UW_Mass_Post.png', dpi=1200)
In [ ]:
In [21]:
figure = triangle.corner(all_mvir, labels=labs, quantiles=[0.16,0.5,0.84], fig=None, weights=Tr.weights,\
plot_contours=True, show_titles=True, title_args={"fontsize": 16}, label_args={"fontsize": 16},\
plot_datapoints=False, bins=20, color='r')
#figure.suptitle("Mass Posterior", fontsize=16, horizontalalignment='left')
In [63]:
figure
Out[63]:
In [62]:
#figure.suptitle("Concentration Posterior", fontsize=16, horizontalalignment='left')
red_patch = mpatches.Patch(color='k')
blue_patch = mpatches.Patch(color='c')
figure.legend(handles=[red_patch, blue_patch], labels=["Cvir from Narrow Mass Bin", "Weighted Triplet Cvir"], fontsize=16)
Out[62]:
In [65]:
figure.savefig('/afs/slac.stanford.edu/u/ki/mwillia1/Thesis/LocalGroupHaloProps/doc/thesis/plots/final_plots/TC_vs_MBC.png', dpi=1200)
In [29]:
h = np.load('/lustre/ki/pfs/mwillia1/LG_project/Consuelo_Boxes/4002/4002hlist.npy')
In [47]:
mw_h = h[np.abs(np.log10(h['mvir'])-11.7) < .1]
m31_h = h[np.abs(np.log10(h['mvir'])-12) < .1]
m33_h = h[np.abs(np.log10(h['mvir'])-10.8) < .1]
In [51]:
mwc = mw_h['rvir']/mw_h['rs']
m31c = m31_h['rvir']/m31_h['rs']
m33c = m33_h['rvir']/m33_h['rs']
mwc.shape, m31c.shape, m33c.shape
Out[51]:
In [52]:
mwc = mwc[0:100000]
m31c = m31c[0:100000]
m33c = m33c[0:100000]
In [58]:
all_mvir = np.transpose(np.vstack((mwc, m31c, m33c)))
labs = ["$Cvir_{MW}$", "$Cvir_{M31}$", "$Cvir_{M33}$"]
In [59]:
figure_c = triangle.corner(all_mvir, labels=labs, quantiles=[0.,0.5,0.84], fig=figure, weights=None,\
plot_contours=True, show_titles=True, title_args={"fontsize": 16}, label_args={"fontsize": 16},\
plot_datapoints=False, bins=20, color='k')
In [ ]:
In [206]:
figure2 = triangle.corner(all_mvir, labels=labs, quantiles=[0.16,0.5,0.84], fig=figure, weights=Tr.weights,\
plot_contours=True, show_titles=True, title_args={"fontsize": 16}, label_args={"fontsize": 16},\
plot_datapoints=False, bins=20, color='r')
In [ ]:
In [209]:
figure
Out[209]:
In [208]:
#figure.suptitle("Mass Posterior", fontsize=16, horizontalalignment='left')
red_patch = mpatches.Patch(color='r')
blue_patch = mpatches.Patch(color='b')
figure.legend(handles=[red_patch, blue_patch], labels=["Triplets", "Pairs"], fontsize=16)
Out[208]:
In [134]:
a=figure.gca()
In [138]:
a.legend_?
In [211]:
figure2.savefig('/afs/slac.stanford.edu/u/ki/mwillia1/Thesis/LocalGroupHaloProps/doc/thesis/plots/final_plots/T_vs_P_mass.png', dpi=1200)
In [48]:
Tr.weights.shape
Out[48]:
In [66]:
xdata = total
plot(xdata, Tr.weights)
xscale('log')
In [67]:
bins = 20
xbins = np.logspace(np.log10(xdata.min()), np.log10(xdata.max()), bins)
ybins = [np.sum(Tr.weights[(xdata[:]>xbins[i]) & (xdata[:]<xbins[i+1])]) for i in range(bins-1)]
print len(xbins)
print len(ybins)
print len(xbins[:bins-1])
In [68]:
plot(xbins[:bins-1], ybins)
xscale('log')
In [40]:
xtrip = np.copy(xbins[:bins-1])
ytrip = np.copy(ybins)
In [45]:
xtripmw = np.copy(xbins[:bins-1])
ytripmw = np.copy(ybins)
In [65]:
xpairmw = np.copy(xbins[:bins-1])
ypairmw = np.copy(ybins)
In [69]:
xpair = np.copy(xbins[:bins-1])
ypair = np.copy(ybins)
In [80]:
plot(xtripmw, ytripmw, color='r')
plot(xpairmw, ypairmw, color='b')
ylabel('P(MMW)')
xlabel('MMW')
red_patch = mpatches.Patch(color='red')
blue_patch = mpatches.Patch(color='blue')
legend(["Trips, M33 dynamics", "Pairs, M33 existence"])
xscale('log')
#savefig("/afs/slac.stanford.edu/u/ki/mwillia1/Thesis/LocalGroupHaloProps/doc/thesis/plots/pdf_mw_plot.png")
In [ ]:
In [ ]:
In [ ]:
In [24]:
all_mvir = np.transpose(np.vstack((np.log10(Tr.MW.Mvir), Tr.MW.Cvir)))
labs = ["$Mvir_{MW}$", "$Cvir_{MW}$"]
In [103]:
figure = triangle.corner(all_mvir, labels=labs, quantiles=[0.16,0.5,0.84], fig=None, weights=Tr.weights,\
plot_contours=True, show_titles=True, title_args={"fontsize": 16}, label_args={"fontsize": 16},\
plot_datapoints=False, bins=20, color='r')
In [54]:
hlist_path = '/lustre/ki/pfs/mwillia1/LG_project/Consuelo_Boxes/4001/4001hlist.npy'
hlist = np.load(hlist_path)
In [55]:
hlist=hlist[np.abs(np.log10(hlist['mvir'])-12)<1.0]
cvir = hlist['rvir']/hlist['rs']
mvir = np.log10(hlist['mvir'])
In [59]:
hlist['mvir']
Out[59]:
In [ ]:
def binMass(mvir, cvir, width):
mvir_bins = np.arange(np.min(mvir), np.max(mvir), width)
In [95]:
bins = np.arange(np.min(mvir), np.max(mvir), .01)
Out[95]:
In [96]:
dat = np.vstack((mvir, cvir)).T
Out[96]:
In [97]:
conc = [[el[1] for el in dat if ((el[0] > bins[i]) & (el[0] < bins[i+1]))] for i in range(len(bins)-1)]
In [98]:
conc_means = np.array([np.mean(np.array(c)) for c in conc])
conc_stds = np.array([np.std(np.array(c)) for c in conc])
In [99]:
conc_means
Out[99]:
In [106]:
faxes = np.reshape(figure.axes, (2,2))
ax=faxes[1,0]
In [107]:
ax.plot(bins[1:], conc_means, lw=2, color='blue')
ax.fill_between(bins[1:], conc_means+conc_stds, conc_means-conc_stds, facecolor='blue', alpha=0.5)
Out[107]:
In [109]:
ax.set_xlim((11,13))
ax.set_ylim((4, 20))
Out[109]:
In [110]:
figure
Out[110]:
In [92]:
hlist_path = '/lustre/ki/pfs/mwillia1/LG_project/Consuelo_Boxes/4001/4001hlist.npy'
hlist = np.load(hlist_path)
hlist=hlist[np.abs(np.log10(hlist['mvir'])-12)<1.0]
In [94]:
scatter(np.log10(hlist['mvir']), hlist['rvir']/hlist['rs'])
Out[94]:
In [111]:
H, xedge, yedge = np.histogram2d(mvir, cvir, bins=[10,10])
In [114]:
H[0,:]
Out[114]:
In [ ]: