<jupyter_start><jupyter_text>Sample 6.2 Multi-parameter Bayesian model: Globular cluster membership<jupyter_code>%matplotlib inline
#%pylab
#posteior distribution of normal data
import numpy as np
import scipy.stats as stats
import matplotlib.pyplot as plt
import astropy.io.fits as fits
import matplotlib
matplotlib.rc('xtick', labelsize=12) 
matplotlib.rc('ytick', labelsize=12) 
from astropy.table import Table
import emcee
import corner

gc = Table.read("N6205.vot",format="votable")

print(gc[0])
ind_gc = (gc['parallax_over_error']>3.) &\
    ( gc['parallax']>-0.03) & ( gc['parallax']<0.2)
print(np.sum(ind_gc))
ind_fld = (gc['parallax_over_error']>3.) &\
    ( gc['parallax']>-0.03) & ( gc['parallax']>=0.2)
print(np.sum(ind_fld))
fig = plt.figure(figsize=(12,5))
ax = fig.add_subplot(121)
ax.plot(gc['ra'][ind_gc],gc['dec'][ind_gc],'k.',markersize=2,alpha=1)
ax.plot(gc['ra'][ind_fld],gc['dec'][ind_fld],'r.',markersize=1,alpha=0.3)
ax.set_xlabel('RA')
ax.set_ylabel('DEC')

ax = fig.add_subplot(122)
ax.plot(gc['pmra'][ind_gc],gc['pmdec'][ind_gc],'k.',markersize=2,alpha=1)
ax.plot(gc['pmra'][ind_fld],gc['pmdec'][ind_fld],'r.',markersize=1,alpha=1)
ax.set_xlim(-14,10)
ax.set_ylim(-14,10)
ax.set_xlabel('pm_RA')
ax.set_ylabel('pm_DE')
ind_gc2 = ind_gc & (gc['pmra']>-4) & (gc['pmra']<-2.5) & (gc['pmdec']>-3) & (gc['pmdec']<-1.7)
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111)
h,xedge = np.histogram(gc['parallax'][ind_gc2 | ind_fld],np.arange(0,5,0.01))
ax.step(xedge[1:],h,where='post')
ax.set_xlim(0,1)
ax.set_xlabel('parallax')
fig = plt.figure(figsize=(10,15))
ax = fig.add_subplot(211)
ax.plot(gc['ra'][ind_gc2],gc['dec'][ind_gc2],'k.',markersize=1)
#ax.plot(gc['ra'][ind_fld],gc['dec'][ind_fld],'r.',markersize=1,alpha=0.3)
ax.set_xlabel('RA')
ax.set_ylabel('DEC')

ax = fig.add_subplot(212)
ax.plot(gc['pmra'][ind_gc2],gc['pmdec'][ind_gc2],'k.',markersize=1)
#ax.plot(gc['pmra'][ind_fld],gc['pmdec'][ind_fld],'r.',markersize=1,alpha=0.3)
# ax.set_xlim(-4,-1)
# ax.set_ylim(-4,-1)
ax.set_xlabel(r'$\mu_{RA}$')
ax.set_ylabel(r'$\mu_{DEC}$')<jupyter_output><empty_output><jupyter_text>estimate the mean proper motions of N6205 $p(\mu,\Sigma|y)\propto p(y|\mu,\Sigma)p(\mu)p(\sigma_\alpha)p(\sigma_\delta)=\frac{\exp(-{\mu}^T\Sigma^{-1}{\mu})}{2\pi\sigma_\alpha^3\sigma_\delta^3\sqrt{1-\rho^2}}$<jupyter_code>y = np.concatenate((gc['pmra'][ind_gc2].data.reshape(np.sum(ind_gc2),1),\
                    gc['pmdec'][ind_gc2].data.reshape(np.sum(ind_gc2),1)),axis=1)
print(y.shape)
print('mean:',np.mean(y,axis=0))
print('std.:',np.std(y,axis=0))
plt.plot(y[:,0],y[:,1],'k.')
def posterior1(theta,y):
    '''
    without 1/sigma**2 prior
    '''
    x = np.zeros_like(y)
    x[:,0] = y[:,0]-theta[0]
    x[:,1] = y[:,1]-theta[1]
    #print(np.mean(y,axis=0),np.mean(x,axis=0))
    s1 = theta[2]
    s2 = theta[3]
    rho = theta[4]
    if np.abs(rho)>1 or s1<=0 or s2<=0  :
        lnp = -1e50
    else:
        lnp = np.sum(-1./(2.*(1-rho**2))*(x[:,0]**2/s1**2-2*rho*x[:,0]*x[:,1]/(s1*s2)+x[:,1]**2/s2**2))-\
            (y.shape[0])*np.log(s1)-(y.shape[0])*np.log(s2)-0.5*y.shape[0]*np.log(1-rho**2)
    return lnp

def posterior2(theta,y):
    '''
    with 1/sigma**2 prior
    '''
    mu = np.zeros_like(y)
    mu[:,0] = y[:,0]-theta[0]
    mu[:,1] = y[:,1]-theta[1]
    s1 = theta[2]
    s2 = theta[3]
    rho = theta[4]
    if np.abs(rho)>1 or s1<=0 or s2<=0:
        lnp = -1e50
    else:
        lnp = np.sum(-1./(2.*(1-rho**2))*(mu[:,0]**2/s1**2-2*rho*mu[:,0]*mu[:,1]/(s1*s2)+mu[:,1]**2/s2**2))-\
            (2.+y.shape[0])*np.log(s1)-(2.+y.shape[0])*np.log(s2)-0.5*y.shape[0]*np.log(1-rho**2)
    return lnp

def mcmc_run(y, p0, lnlikeli = posterior1, ndim = 2, nwalkers = 50):
    

    sampler = emcee.EnsembleSampler(nwalkers, \
        ndim, lnlikeli, \
        args=[y])

    pos, prob, state = sampler.run_mcmc(p0, 100)
    sampler.reset()

    sampler.run_mcmc(pos, 1000)

    samples = sampler.chain[:, :, :].reshape((-1, ndim))
    return samples

    
## always test your ln posterior before MCMC
theta1= [-2.5,-2.5,0.2,0.3,0.]
print(posterior1(theta1,y),posterior2(theta1,y))
theta2= [-0.7,-0.1,0.2,0.3,0.]
print(posterior1(theta2,y),posterior2(theta2,y))
## using posterior1 without 1/sigma**2 prior

nwalkers = 50
ndim = 5

p0=np.zeros((nwalkers,ndim))
p0[:,0] = np.random.rand(nwalkers)-3.
p0[:,1] = np.random.rand(nwalkers)-3.
p0[:,2] = np.random.rand(nwalkers)*0.3
p0[:,3] = np.random.rand(nwalkers)*0.3
p0[:,4] = np.random.rand(nwalkers)*2.-1.

samples1 = mcmc_run(y, p0, lnlikeli = posterior1, ndim = ndim, nwalkers = nwalkers)

fig = corner.corner(samples1,\
            labels=[r"$\mu_\alpha$",r"$\mu_\delta$", r"$\sigma_\alpha$", r"$\sigma_\delta$",r"$\rho$"],\
            quantiles=[0.16, 0.5, 0.84],\
            show_titles=True, \
            title_kwargs={"fontsize": 12})
## using posterior2 with 1/sigma**2 prior

samples2 = mcmc_run(y, p0, lnlikeli = posterior2, ndim = ndim, nwalkers = nwalkers)

fig = corner.corner(samples2,\
            labels=[r"$\mu_\alpha$",r"$\mu_\delta$", r"$\sigma_\alpha$", r"$\sigma_\delta$",r"$\rho$"],\
            quantiles=[0.16, 0.5, 0.84],\
            show_titles=True, \
            title_kwargs={"fontsize": 12})<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.24737945492662475, 'fraction_numerical': 0.06041547550981513, 'mean_word_length': 4.355102040816327, 'pattern_counts': {'":': 2, '<': 16, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '2', 'dir': 'jupyter-structured-clean-dedup', 'id': '493492', 'license': 'MIT', 'n_tokens_mistral': 2503, 'n_tokens_neox': 2443, 'n_words': 289, 'path': 'sample6.2_MultiParamModel_GCpropermotion.ipynb', 'repo_name': 'liuchaonaoc/Astrostatistics2022', 'size': '5247'}