<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Mass Maps From Mass-Luminosity Inference Posterior
In this notebook we start to explore the potential of using a mass-luminosity relation posterior to refine mass maps. 
Content
Step2:  Probability Functions
Step3:  Results
Step4:  Turning into Probabilistic Catalogue
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import rc
rc('text', usetex=True)
from bigmali.grid import Grid
from bigmali.prior import TinkerPrior
from bigmali.hyperparameter import get
import numpy as np
from scipy.stats import lognorm
from numpy.random import normal
#globals that functions rely on
grid = Grid()
prior = TinkerPrior(grid)
a_seed = get()[:-1]
S_seed = get()[-1]
mass_points = prior.fetch(grid.snap(0)).mass[2:-2] # cut edges
tmp = np.loadtxt('/Users/user/Code/PanglossNotebooks/MassLuminosityProject/SummerResearch/mass_mapping.txt')
z_data = tmp[:,0]
lobs_data = tmp[:,1]
mass_data = tmp[:,2]
ra_data = tmp[:,3]
dec_data = tmp[:,4]
sigobs = 0.05
def fast_lognormal(mu, sigma, x):
    return  (1/(x * sigma * np.sqrt(2 * np.pi))) * np.exp(- 0.5 * (np.log(x) - np.log(mu)) ** 2 / sigma ** 2)
Explanation: Mass Maps From Mass-Luminosity Inference Posterior
In this notebook we start to explore the potential of using a mass-luminosity relation posterior to refine mass maps. 
Content:
- Math
- Imports, Constants, Utils, Data
- Probability Functions
- Results
- Discussion
Math
Infering mass from mass-luminosity relation posterior ...
\begin{align}
P(M|L_{obs},z,\sigma_L^{obs}) &= \iint P(M|\alpha, S, L_{obs}, z)P(\alpha, S|L_{obs},z,\sigma_L^{obs})\ d\alpha dS\
&\propto \iiint P(L_{obs}| L,\sigma_L^{obs})P(L|M,\alpha,S,z)P(M|z)P(\alpha, S|L_{obs},z,\sigma_L^{obs})\ dLd\alpha dS\
&\approx \frac{P(M|z)}{n_{\alpha,S}}\sum_{\alpha,S \sim P(\alpha, S|L_{obs},z,\sigma_L^{obs})}\left( \frac{1}{n_L}\sum_{L\sim P(L|M,\alpha,S,z)}P(L_{obs}|L,\sigma_L^{obs})\right)\
&= \frac{P(M|z)}{n_{\alpha,S}}\sum_{\alpha,S \sim P(\alpha, S|L_{obs},z,\sigma_L^{obs})}f(M;\alpha,S,z)\
\end{align}
Refine for individual halo ...
\begin{align}
P(M_k|L_{obs},z,\sigma_L^{obs}) &= \iint P(M_k|\alpha, S, L_{obs\ k}, z_k)P(\alpha, S|L_{obs},z,\sigma_L^{obs})\ d\alpha dS\
&\propto \iiint P(L_{obs\ k}| L_k,\sigma_L^{obs})P(L_k|M_k,\alpha,S,z_k)P(M_k|z_k)P(\alpha, S|L_{obs},z,\sigma_L^{obs})\ dLd\alpha dS\
&\approx \frac{P(M_k|z_k)}{n_{\alpha,S}}\sum_{\alpha,S \sim P(\alpha, S|L_{obs},z,\sigma_L^{obs})}\left( \frac{1}{n_L}\sum_{L\sim P(L_k|M_k,\alpha,S,z_k)}P(L_{obs\ k}|L_k,\sigma_L^{obs})\right)\
&=\frac{P(M_k|z_k)}{n_{\alpha,S}}\sum_{\alpha,S \sim P(\alpha, S|L_{obs},z,\sigma_L^{obs})}f(M_k;\alpha,S,z_k)\
\end{align}
Can also factor it more conventionally for MCMC ...
\begin{align}
\underbrace{P(M_k|L_{obs},z,\sigma_L^{obs})}{posterior} 
&\propto \underbrace{P(M_k|z_k)}{prior}\underbrace{\iiint P(L_{obs\ k}| L_k,\sigma_L^{obs})P(L_k|M_k,\alpha,S,z_k)P(\alpha, S|L_{obs},z,\sigma_L^{obs})\ dLd\alpha dS}_{likelihood}\
\end{align}
In the code we have the following naming convention:
- p1 for $P(M|z)$
- p2 for $P(\alpha, S|L_{obs},z,\sigma_L^{obs})$
- p3 for $P(L_k|M_k,\alpha,S,z_k)$
- p4 for $P(L_{obs\ k}|L_k, \sigma^{obs}_L)$
We use the terms eval and samp to help distinguish between evaluating a distribution and sampling from it. 
Imports, Constants, Utils, Data
End of explanation
def p1_eval(zk):
    return prior.fetch(grid.snap(zk)).prob[2:-2]
def p2_samp(nas=100):
    
    a is fixed on hyperseed,
    S is normal distribution centered at hyperseed.
    
    return normal(S_seed, S_seed / 10, size=nas)
def p3_samp(mk, a, S, zk, nl=100):
    mu_lum = np.exp(a[0]) * ((mk / a[2]) ** a[1]) * ((1 + zk) ** (a[3]))
    return lognorm(S, scale=mu_lum).rvs(nl)
    
def p4_eval(lobsk, lk, sigobs):
    return fast_lognormal(lk, sigobs, lobsk)
    
def f(a, S, zk, lobsk, nl=100):
    ans = []
    for mk in mass_points:
        tot = 0
        for x in p3_samp(mk, a, S, zk, nl):
            tot += p4_eval(lobsk, x, sigobs)
        ans.append(tot / nl)
    return ans
def mass_dist(ind=1, nas=10, nl=100):
    lobsk = lobs_data[ind]
    zk = z_data[ind]
    tot = np.zeros(len(mass_points))
    for S in p2_samp(nas):
        tot += f(a_seed, S, zk, lobsk, nl)
    prop = p1_eval(zk) * tot / nas
    return prop / np.trapz(prop, x=mass_points)
Explanation: Probability Functions
End of explanation
plt.subplot(3,3,1)
dist = p1_eval(zk)
plt.plot(mass_points, dist)
plt.gca().set_xscale('log')
plt.gca().set_yscale('log')
plt.ylim([10**-25, 10])
plt.xlim([mass_points.min(), mass_points.max()])
plt.title('Prior')
plt.xlabel(r'Mass $(M_\odot)$')
plt.ylabel('Density')
for ind in range(2,9):
    plt.subplot(3,3,ind)
    dist = mass_dist(ind)
    plt.plot(mass_points, dist, alpha=0.6, linewidth=2)
    plt.xlim([mass_points.min(), mass_points.max()])
    plt.gca().set_xscale('log')
    plt.gca().set_yscale('log')
    plt.ylim([10**-25, 10])
    plt.gca().axvline(mass_data[ind], color='red', linewidth=2, alpha=0.6)
    plt.title('Mass Distribution')
    plt.xlabel(r'Mass $(M_\odot)$')
    plt.ylabel('Density')
    
    
# most massive
ind = np.argmax(mass_data)
plt.subplot(3,3,9)
dist = mass_dist(ind)
plt.plot(mass_points, dist, alpha=0.6, linewidth=2)
plt.gca().set_xscale('log')
plt.gca().set_yscale('log')
plt.xlim([mass_points.min(), mass_points.max()])
plt.ylim([10**-25, 10])
plt.gca().axvline(mass_data[ind], color='red', linewidth=2, alpha=0.6)
plt.title('Mass Distribution')
plt.xlabel(r'Mass $(M_\odot)$')
plt.ylabel('Density')
# plt.tight_layout()
plt.gcf().set_size_inches((10,6))
Explanation: Results
End of explanation
index = range(2,9) + [np.argmax(mass_data)]
plt.title('Simple Sketch of Field of View')
plt.scatter(ra_data[index], dec_data[index] , s=np.log(mass_data[index]), alpha=0.6)
plt.xlabel('ra')
plt.ylabel('dec');
Explanation: Turning into Probabilistic Catalogue
End of explanation
<END_TASK>