<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Inferring Cluster Model Parameters from an X-ray Image
Forward modeling is always instructive
Step1:  This PGM illustrates the joint PDF for the parameters and the data, which can be factorised as
Step2:  Good. Here's the code that is being run, inside the "XrayData" class
Step3:  To normalize this, we need to take care not to try and exponentiate any very large or small numbers...
Step4:  Let's plot this as a 2D probability density map.
<ASSISTANT_TASK:>
Python Code:
# import cluster_pgm
# cluster_pgm.inverse()
from IPython.display import Image
Image(filename="cluster_pgm_inverse.png")
Explanation: Inferring Cluster Model Parameters from an X-ray Image
Forward modeling is always instructive: we got a good sense of the parameters of our cluster + background model simply by generating mock data and visualizing it.
The "inverse problem", also known as "inference," is to learn the parameters of an assumed model from a set of data. Intuitively we can see how it is going to work: try a lot of possible parameter combinations, and see which ones "match" the data.
Our inability to guess parameter values accurately first time shows that we are uncertain about them. In Bayesian inference, we use probability distributions to describe this uncertainty mathematically.
The sampling distribution ${\rm Pr}(d|\theta,H)$ encodes uncertainty about what might have been, given a model (or hypothesis) $H$ with parameters $\theta$. It allows us to generate mock datasets that are similar to the data that we do observe.
Before we take any data, our uncertainty about our model parameter values is encoded in the prior PDF for the parameters given the model, ${\rm Pr}(\theta|H)$.
Similarly, the sampling distribution ${\rm Pr}(d|\theta,H)$ can be thought of as the prior over datasets - a PDF could be imagined and used to generate mock data without us ever having seen any real data at all!
Probability
The idea of using probability distributions to quantify the uncertainty in our model parameters (and indeed in the models themselves) is due to Pierre Simon Laplace (1774), who rediscovered Thomas Bayes' earlier results on the probability of future events given their past history.
Let's remind ourselves how probabilities work
Laplace and Bayes' key result is the following, usually referred to as "Bayes' Theorem:"
${\rm Pr}(\theta|d,H) = \frac{1}{{\rm Pr}(d|H)}\;{\rm Pr}(d|\theta,H)\;{\rm Pr}(\theta|H)$
What you know about your model parameters given the data is what you knew about them before $\left[ {\rm Pr}(\theta|H) \right]$, combined with what the data are telling you $\left[ {\rm Pr}(d|\theta,H) \right]$.
${\rm Pr}(\theta|d,H)$ is called the posterior probability distribution for the parameters given the data, and is the general solution to the inverse problem.
Both the posterior and prior PDFs are functions of the model parameters. The sampling distribution is a function of the data given the parameters - when written as a function of $\theta$ it is called the likelihood of the parameters given the model.
PGMs for Inverse Problems
Here's the probabilistic graphical model for the inverse X-ray cluster model problem.
Q: Spot the difference!
End of explanation
%load_ext autoreload
%autoreload 2
import cluster
lets = cluster.XrayData()
lets.read_in_data()
lets.set_up_maps()
x0,y0 = 328,328    # The center of the image is 328,328
S0,b = 0.001,1e-6  # Cluster and background surface brightness, arbitrary units
beta = 2.0/3.0     # Canonical value is beta = 2/3
rc = 12            # Core radius, in pixels
logprob = lets.evaluate_unnormalised_log_posterior(x0,y0,S0,rc,beta,b)
print logprob
Explanation: This PGM illustrates the joint PDF for the parameters and the data, which can be factorised as:
$\prod_k \; {\rm Pr}(N_k\;|\;\mu_k(\theta),{\rm ex}_k,{\rm pb}_k,H) \; {\rm Pr}(\,\theta\,|H)$
It can also be factorised to:
${\rm Pr}(\,\theta\,|{N_k}\,H) \; {\rm Pr}({N_k}\,|H)$
which is, up to the normalizing constant, the posterior PDF for the model parameters, given all the data ${N_k}$.
PGMs can be used to design inferences
Calculating Posterior PDFs
Notice that the prior PDF ${\rm Pr}(\theta|H)$ and the likelihood function ${\rm Pr}(d|\theta,H)$ can typically be evaluated at any point in the parameter space.
This means that we can always simply evaluate the posterior PDF on a grid (or at least attempt to), and normalize it by numerical integration.
Let's do this for a simplified version of our X-ray cluster model.
End of explanation
import numpy as np
npix = 15
xmin,xmax = 327.7,328.3
ymin,ymax = 346.4,347.0
x0grid = np.linspace(xmin,xmax,npix)
y0grid = np.linspace(ymin,ymax,npix)
logprob = np.zeros([npix,npix])
for i,x0 in enumerate(x0grid):
    for j,y0 in enumerate(y0grid):
        logprob[j,i] = lets.evaluate_unnormalised_log_posterior(x0,y0,S0,rc,beta,b)
    print "Done column",i
print logprob[0:5,0]
Explanation: Good. Here's the code that is being run, inside the "XrayData" class:
```python
    def evaluate_log_prior(self):
        # Uniform in all parameters...
        return 0.0
def evaluate_log_likelihood(self):
    self.make_mean_image()
    # Return un-normalized Poisson sampling distribution:
    # log (\mu^N e^{-\mu} / N!) = N log \mu - \mu + constant
    return np.sum(self.im * np.log(self.mu) - self.mu)
def evaluate_unnormalised_log_posterior(self,x0,y0,S0,rc,beta,b):
    self.set_pars(x0,y0,S0,rc,beta,b)
    return self.evaluate_log_likelihood() + self.evaluate_log_prior()
```
Now let's try evaluating the 2D posterior PDF for cluster position, conditioned on reasonable values of the cluster and background flux, cluster size and beta:
End of explanation
Z = np.max(logprob)
prob = np.exp(logprob - Z)
norm = np.sum(prob)
prob /= norm
print prob[0:5,0]
Explanation: To normalize this, we need to take care not to try and exponentiate any very large or small numbers...
End of explanation
import astropy.visualization as viz
import matplotlib.pyplot as plt
%matplotlib inline
plt.rcParams['figure.figsize'] = (10.0, 10.0)
plt.imshow(prob, origin='lower', cmap='Blues', interpolation='gaussian', extent=[xmin,xmax,ymin,ymax])
plt.xlabel('x / pixels')
plt.ylabel('y / pixels')
Explanation: Let's plot this as a 2D probability density map.
End of explanation
<END_TASK>