In [1]:
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import linmix
from astropy.table import Table
import astropy.io.ascii as ascii
import corner
from lifelines import KaplanMeierFitter
from matplotlib import rc
from matplotlib.ticker import MultipleLocator
rc('font', **{'family': 'serif', 'serif': ['Computer Modern']})
rc('text', usetex=True)
#%matplotlib inline
%matplotlib qt5
In [5]:
Tab=Table.read("Table.fit") #Taurus fluxes and disk mass Tab
Tab_Mass=Table.read("Table_Mass.fit") #Taurus stellar masses (2-3 models) Tab
Lup_Tab=Table.read('Lupus_Tab.fit') #Lupus Tab
ChaI_Tab=Table.read('ChaI_Tab.fit')
UppSco_Tab=Table.read('UppSco_Tab.fit')
Md_Lup=Lup_Tab['Md']
e_Md_Lup=Lup_Tab['err_Md']
Md_Lup[notd_Lup]=3*e_Md_Lup[notd_Lup] #Value for upperlimits
In [47]:
Tau_20 = ascii.read('Md_20_vs_Ms.pyout')
Tau2_20 = ascii.read('Md_20_vs_Ms2.pyout')
Tau3_20 = ascii.read('Md_20_vs_Ms3.pyout')
Lup_20 = ascii.read('Lup_Md_vs_Ms.pyout')
ChaI_20 = ascii.read('ChaI_Md_vs_Ms.pyout')
UppSco_20 = ascii.read('UppSco_Md_vs_Ms.pyout')
a_Tau=Tau_20['alpha']
b_Tau=Tau_20['beta']
sig_Tau=np.sqrt(Tau_20['sigsqr'])
corr_Tau=np.sqrt(Tau_20['corr'])
a_Lup=Lup_20['alpha']
b_Lup=Lup_20['beta']
sig_Lup=np.sqrt(Lup_20['sigsqr'])
corr_Lup=np.sqrt(Lup_20['corr'])
a_ChaI=ChaI_20['alpha']
b_ChaI=ChaI_20['beta']
sig_ChaI=np.sqrt(ChaI_20['sigsqr'])
corr_ChaI=np.sqrt(ChaI_20['corr'])
a_UppSco=UppSco_20['alpha']
b_UppSco=UppSco_20['beta']
sig_UppSco=np.sqrt(UppSco_20['sigsqr'])
corr_UppSco=np.sqrt(UppSco_20['corr'])
d_Tau=Tab['l']!='<' #observed sources
d_Lup=Lup_Tab['l_Md']!='<' #observed sources
d_ChaI=ChaI_Tab['l_Md']!='<' #observed sources
d_UppSco=UppSco_Tab['l_Md']!='<' #observed sources
notd_Lup=np.logical_not(d_Lup)
logMs_Tau=Tab['LogM*']
logMd20_Tau=Tab['LogMd_20']
logMs_Lup=np.log10(Lup_Tab['M*'])
logMd_Lup=np.log10(Lup_Tab['Md'])
logMs_ChaI=ChaI_Tab['logM*']
logMd_ChaI=np.log10(ChaI_Tab['Md'])
logMs_UppSco=UppSco_Tab['logM*']
logMd_UppSco=np.log10(UppSco_Tab['Md_20'])
delta=(d_Tau,d_Lup,d_ChaI,d_UppSco)
#notdelta = np.logical_not(delta) #upperlimits
notdelta=(np.logical_not(d_Tau),np.logical_not(d_Lup),np.logical_not(d_ChaI),np.logical_not(d_UppSco))
a_tot=(a_Tau,a_Lup,a_ChaI,a_UppSco)
b_tot=(b_Tau,b_Lup,b_ChaI,b_UppSco)
s_tot=(sig_Tau,sig_Lup,sig_ChaI,sig_UppSco)
corr_tot=(corr_Tau,corr_Lup,corr_ChaI,corr_UppSco)
logMs=(logMs_Tau,logMs_Lup,logMs_ChaI,logMs_UppSco)
logMd=(logMd20_Tau,logMd_Lup,logMd_ChaI,logMd_UppSco)
a_mean=(a_Tau.mean(),a_Lup.mean(),a_ChaI.mean(),a_UppSco.mean())
a_std=(a_Tau.std(),a_Lup.std(),a_ChaI.std(),a_UppSco.std())
b_mean=(b_Tau.mean(),b_Lup.mean(),b_ChaI.mean(),b_UppSco.mean())
b_std=(b_Tau.std(),b_Lup.std(),b_ChaI.std(),b_UppSco.std())
sig_mean=(sig_Tau.mean(),sig_Lup.mean(),sig_ChaI.mean(),sig_UppSco.mean())
sig_std=(sig_Tau.std(),sig_Lup.std(),sig_ChaI.std(),sig_UppSco.std())
print a_mean
print a_std
print b_mean
print sig_mean
In [13]:
fig, ax = plt.subplots(2, 3, figsize=(5, 5) ,sharey=True,sharex=True)
plt.subplots_adjust(wspace=1e-20,hspace=1e-20)
survey=('Taurus','Lupus','Cham I','Upper Sco')
color=('magenta','blue','orange','green')
ax[1,0].set_xlabel('log(M$_{_*}$/M$_\odot$)', fontsize='16')
ax[1,1].set_xlabel('log(M$_{_*}$/M$_\odot$)', fontsize='16')
ax[1,2].set_xlabel('log(M$_{_*}$/M$_\odot$)', fontsize='16')
#ax[0,3].set_xlabel('log(M$_{_*}$/M$_\odot$)', fontsize='16')
ax[0,0].set_ylabel('log(M$_d$/M$_\oplus$)', fontsize='16')
ax[1,0].set_ylabel('log(M$_d$/M$_\oplus$)', fontsize='16')
for j in range(3):
ax[0,j].text(0, -1.5, survey[j],color=color[j], fontsize=15)
ax[0,j].text(-1.3, 2.6, r'$<\alpha>$='"{0:0.2f}".format(a_mean[j]),color=color[j], fontsize=13)
ax[0,j].text(-1.3, 2.3 , r'$<\beta>$='"{0:0.2f}".format(b_mean[j]),color=color[j], fontsize=13)
ax[0,j].text(-1.3, 2, r'$<\sigma>$='"{0:0.2f}".format(sig_mean[j]),color=color[j], fontsize=13)
ax[0,j].set_ylim(-2,2.9)
ax[0,j].set_xlim(-1.4,0.6)
ax[0,j].errorbar(logMs[j][delta[j]],logMd[j][delta[j]] ,fmt='o', color=color[j], alpha=0.5)
ax[0,j].errorbar(logMs[j][notdelta[j]],logMd[j][notdelta[j]],fmt='v',color='grey', alpha=0.3)
ax[0,j].tick_params(axis='x', labelsize=12)
ax[0,j].tick_params(axis='y', labelsize=12)
ax[1,0].text(0, -1.5, survey[3],color=color1[3], fontsize=15)
ax[1,0].text(-1.3, 2.6, r'$<\alpha>$='"{0:0.2f}".format(a_mean[3]),color=color[3], fontsize=13)
ax[1,0].text(-1.3, 2.3 , r'$<\beta>$='"{0:0.2f}".format(b_mean[3]),color=color[3], fontsize=13)
ax[1,0].text(-1.3, 2, r'$<\sigma>$='"{0:0.2f}".format(sig_mean[3]),color=color[3], fontsize=13)
ax[1,0].set_ylim(-2,2.9)
ax[1,0].set_xlim(-1.4,0.6)
ax[1,0].errorbar(logMs[3][delta[3]],logMd[3][delta[3]] ,fmt='o', color=color[3], alpha=0.5)
ax[1,0].errorbar(logMs[3][notdelta[3]],logMd[3][notdelta[3]],fmt='v',color='grey', alpha=0.3)
ax[1,0].tick_params(axis='x', labelsize=12)
ax[1,0].tick_params(axis='y', labelsize=12)
for i in xrange(0, len(a_tot[j]), 1000):
xs = np.arange(-10,10)
ys1 = a_tot[j][i] + xs * b_tot[j][i]
ax[0,j].plot(xs,ys1, color=color1[j], alpha=0.02)
ys3= a_tot[3][i] + xs * b_tot[3][i]
ax[1,0].plot(xs,ys3, color=color1[3], alpha=0.02)
ax[1,2].plot(xs,a_mean[j]+xs*b_mean[j], color=color[j], alpha=0.02)
ax[1,2].plot(xs,a_mean[3]+xs*b_mean[3], color=color[3], alpha=0.02)
plt.show()
Histograms
In [63]:
w1 = np.ones_like(a_tot[0])/float(len(a_tot[0]))
fig1, ax1 = plt.subplots(1, 4, figsize=(5, 5) ,sharey=True)
plt.subplots_adjust(wspace=0,hspace=0.01)
label=(r'$\alpha$',r'$\beta$','RMS','Corr')
hatch=("\\\ ","//","\\ ",'///')
ax1[0].set_ylabel('P(X|M$_{_*}$,M$_d$)', fontsize=16)
ax1[0].set_xlim(0.2,2.4)
ax1[1].set_xlim(0.7,3.8)
ax1[2].set_xlim(0.3,1.05)
ax1[3].set_xlim(-0.1,1.05)
bin1=np.arange(0.2,2.4,0.05)
bin2=np.arange(0.5,3.8,0.05)
bin3=np.arange(0.3,1.05,0.05)
bin4=np.arange(-0.1,1.05,0.05)
for j in range(4):
ax1[j].tick_params(axis='x', labelsize=12)
ax1[j].tick_params(axis='y', labelsize=12)
ax1[j].set_ylim(0,0.55)
ax1[j].set_xlabel(label[j],fontsize=14)
ax1[0].hist(a_tot[j],bins=bin1 ,weights=w1, edgecolor=color[j], hatch=hatch[j] ,histtype='step',label=survey[j])
ax1[1].hist(b_tot[j],bins=bin2 ,weights=w1,edgecolor=color[j], hatch=hatch[j] ,histtype='step',label=survey[j])
ax1[2].hist(s_tot[j],bins=bin3 ,weights=w1,edgecolor=color[j], hatch=hatch[j] ,histtype='step',label=survey[j])
ax1[3].hist(corr_tot[j],bins=bin4 ,weights=w1,edgecolor=color[j], hatch=hatch[j],histtype='step',label=survey[j])
ax1[0].legend(loc='upper left',prop={'size':10})
plt.show()
Survival function with Kaplan-Meier estimator
In [31]:
plot='ratios'
#plot='diskmass'
if plot=='ratios': #survival distrib for ratios Mdisk/Mstellar_host
M_Tau = (10**Tab['LogMd_20'])/(10**Tab['LogM*'])
M_Tau2 = (10**Tab['LogMd_20'])/(10**Tab_Mass['LogM*2'])
M_Tau3 = (10**Tab['LogMd_20'])/(10**Tab_Mass['LogM*3'])
M_Lup = Lup_Tab['Md']/(Lup_Tab['M*'])
M_ChaI = ChaI_Tab['Md']/(10**ChaI_Tab['logM*'])
M_UppSco = UppSco_Tab['Md_20']/(10**UppSco_Tab['logM*'])
if plot=='diskmass': #survival distrib for Disk mass
M_Tau = 10**Tab['LogMd_20']
M_Tau2 = 10**Tab['LogMd_20']
M_Tau3 = 10**Tab['LogMd_20']
M_Lup = Lup_Tab['Md']*1.
M_ChaI = ChaI_Tab['Md']*1.
M_UppSco = UppSco_Tab['Md_20']*1.
C=delta
kmf_Tau= KaplanMeierFitter()
kmf_Lup= KaplanMeierFitter()
kmf_ChaI= KaplanMeierFitter()
kmf_UppSco= KaplanMeierFitter()
kmf_Tau.fit(M_Tau, event_observed=C[0],left_censorship=True,label="kmf",alpha=0.5)
kmf_Lup.fit(M_Lup, event_observed=C[1],left_censorship=True,label="kmf",alpha=0.5)
kmf_ChaI.fit(M_ChaI, event_observed=C[2],left_censorship=True,label="kmf",alpha=0.5)
kmf_UppSco.fit(M_UppSco, event_observed=C[3],left_censorship=True,label="kmf",alpha=0.5)
s_Tau=1-kmf_Tau.cumulative_density_.values.ravel()
s_Lup=1-kmf_Lup.cumulative_density_.values.ravel()
s_ChaI=1-kmf_ChaI.cumulative_density_.values.ravel()
s_UppSco=1-kmf_UppSco.cumulative_density_.values.ravel()
t_Tau=kmf_Tau.timeline
t_Lup=kmf_Lup.timeline
t_ChaI=kmf_ChaI.timeline
t_UppSco=kmf_UppSco.timeline
Dp_Tau=kmf_Tau.confidence_interval_['kmf_upper_0.50']-(1-s_Tau)
Dn_Tau=(1-s_Tau)-kmf_Tau.confidence_interval_['kmf_lower_0.50']
Dp_Lup=kmf_Lup.confidence_interval_['kmf_upper_0.50']-(1-s_Lup)
Dn_Lup=(1-s_Lup)-kmf_Lup.confidence_interval_['kmf_lower_0.50']
Dp_ChaI=kmf_ChaI.confidence_interval_['kmf_upper_0.50']-(1-s_ChaI)
Dn_ChaI=(1-s_ChaI)-kmf_ChaI.confidence_interval_['kmf_lower_0.50']
Dp_UppSco=kmf_UppSco.confidence_interval_['kmf_upper_0.50']-(1-s_UppSco)
Dn_UppSco=(1-s_UppSco)-kmf_UppSco.confidence_interval_['kmf_lower_0.50']
t=(t_Tau,t_Lup,t_ChaI,t_UppSco)
s=(s_Tau,s_Lup,s_ChaI,s_UppSco)
Dn=(Dn_Tau,Dn_Lup,Dn_ChaI,Dn_UppSco)
Dp=(Dp_Tau,Dp_Lup,Dp_ChaI,Dp_UppSco)
ax = plt.subplot(111)
for i in range(4):
ax.errorbar(t[i],s[i],yerr=[Dn[i],Dp[i]],drawstyle='steps-post',linewidth=3,
color=color[i],label=survey[i],elinewidth=1)
minorLocator = MultipleLocator(0.05)
ax.yaxis.set_minor_locator(minorLocator)
plt.legend()
plt.xscale('log')
plt.ylim(0, 1)
if plot=='ratios':
plt.ylabel('f($>$M$_d$/M$_{_*}$)', fontsize=16)
plt.xlabel('M$_d$/M$_{_*}$',fontsize=16)
if plot=='diskmass':
plt.ylabel('f($>$M$_d$)', fontsize=16)
plt.xlabel('M$_d$',fontsize=16)
plt.show()