# Bayesian Model Fitting in Python with emcee
## *Computational Research Day 2018 Workshop*


***
By AM Geller and AA Miller <br/>
April 10, 2018

Today we will focus on a seemingly simple problem: fitting a straight line to data.

Though we've all done this before, we've probably at some point oversimplified the analysis or at least not fully understood all the assumptions that went into our analysis. 

This workshop is inspired by the paper [Data Analysis Recipes: Fitting a Model to Data](https://arxiv.org/abs/1008.4686) by Hogg, Bovy, & Lang. If you haven't already read this paper, I ***highly*** recommend it.


```python
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Ellipse
from scipy.optimize import minimize

import emcee
import corner

#%matplotlib notebook
%matplotlib inline

np.random.seed(seed=2222)
ncores = 6 # adjust this to match your machine
```

The first thing we need is the data.  Let's start with a prepared data set from that paper.  Later we can move to a real data set (maybe your own).

In this data set, you'll notice that each $\left(x,y\right)$ point has an uncertainty characterized by a guassian with variance $\left(\sigma_x^, \sigma_y^2\right)$, while there is also covariance between the $x$ and $y$ measurements described by a correlation coefficient $\rho_{xy}$.  I found [this nice 1-page info sheet online](http://www.cs.utah.edu/~tch/CS4300/resources/refs/ErrorEllipses.pdf) if you want a further reference.

The full uncertainty covariance matrix for each data point is given by:

$$S_i = \left[ {\begin{array}{cc}
        \sigma_{xi}^2 & \rho_{xyi}\sigma_{xi}\sigma_{yi} \\
        \rho_{xyi}\sigma_{xi}\sigma_{yi} & \sigma_{yi}^2  \\
        \end{array}
       }
       \right] 
       $$


```python
obs_x = np.array([203., 58., 210., 202., 198., 158., 
                165., 201., 157., 131., 166., 160., 
                186., 125., 218., 146.])
obs_y = np.array([495., 173., 479., 504., 510., 416., 
                393., 442., 317., 311., 400., 337., 
                423., 334., 533., 344.])
sigma_x = np.array([5., 9., 4., 4., 11., 7., 
                    5., 5., 5., 6., 6., 5., 
                    9., 8., 6., 5.])
sigma_y = np.array([21., 15., 27., 14., 30., 16., 
                    14., 25., 52., 16., 34., 31., 
                    42., 26., 16., 22.])
rho_xy = np.array([-0.33, 0.67, -0.02, -0.05, -0.84, -0.69,
                    0.30, -0.46, -0.03, 0.50, 0.73, -0.52, 
                    0.90, 0.40, -0.78, -0.56])
```

Let's plot the data.


```python
def plotData(ax, x, y, xerr, yerr, rhoxy, zorder = 10):
    ax.scatter(x, y, color = 'black', zorder = 10)
    for (xx, yy, sx, sy, rxy) in zip(x,y,xerr, yerr, rhoxy):
        cov = np.array([[ sx**2, rxy*sx*sy],
                        [rxy*sx*sy, sy**2.]])
        w, v  = np.linalg.eig(cov)
        theta = np.arctan(2. * rxy * sx * sy / (sx**2. - sy**2.))/2.
        ax.add_artist(Ellipse((xx, yy), 2.*w[0]**0.5, 2.*w[1]**0.5, \
                              angle = np.rad2deg(theta),\
                              facecolor="none", edgecolor="black", zorder = zorder))
        
f, ax = plt.subplots(1, figsize = (5,5))
plotData(ax, obs_x, obs_y, sigma_x, sigma_y, rho_xy)
ax.set_xlim(0,300)
ax.set_ylim(0,700)
ax.set_xlabel('x')
ax.set_ylabel('y')
```

Now, let's fit the line :

$y = mx + b$.  

How would you do this?

Probably the first attempt would be some kind of least squares optimization ... but ...

A standard assumption of the least-squares method is that the independent variable (typically $x$) is measured with perfect precision. Thus, standard least squares cannot account for the uncertainties on both $x$ and $y$ (simultaneously).  As a demonstration of ***WHAT NOT TO DO*** (but is sometimes attempted), here's an example of a "forward-reverse" fitting procedure, where you first fit by ignoring $\sigma_x$ then do the inverse and ignore $\sigma_y$.  

I repeat ***DO NOT DO THIS*** for a real publication.  This is only to illustrate a point.


```python
f, ax = plt.subplots(1, figsize = (5,5))
plotData(ax, obs_x, obs_y, sigma_x, sigma_y, rho_xy)
ax.set_xlim(0,300)
ax.set_ylim(0,700)
ax.set_xlabel('x')
ax.set_ylabel('y')

#"abominable" forward–reverse fitting procedure
#DON'T DO THIS!
#polyfit with uncertainties on y ("forward")
p_xy, pcov = np.polyfit(obs_x, obs_y, 1, w = 1./sigma_y, cov = True)
err_xy = np.sqrt(np.diag(pcov))
p_xy_eval = np.poly1d(p_xy)
print("The best-fit value for the ('forward') slope and intercept are: {:.4f} +/- {:.4f} and {:.4f} +/- {:.4f}"\
      .format(p_xy[0], err_xy[0], p_xy[1], err_xy[1]))
ax.plot([0,300],p_xy_eval([0,300]),color='C0')

#polyfit with uncertainties on x ("reverse")
p_yx, pcov = np.polyfit(obs_y, obs_x, 1, w = 1./sigma_x, cov = True)
err_yx = np.sqrt(np.diag(pcov))
p_yx_eval = np.poly1d(p_yx)
inv = [1./p_yx[0], -p_yx[1]/p_yx[0]]
inv_err = [inv[0]**2.*err_yx[0],  (p_yx[0]**-2 * err_yx[1]**2. + (p_yx[1]/p_yx[0]**2.)**2. * err_yx[0]**2.)**0.5]
print("The best-fit value for the ('reverse') slope and intercept are: {:.4f} +/- {:.4f} and {:.4f} +/- {:.4f}"\
      .format(inv[0], inv_err[0], inv[1], inv_err[1]))
ax.plot(p_yx_eval([0.,700.]), [0.,700.],'--',color='C1')

mFreq = p_xy[0]
bFreq = p_xy[1]

```

Again ***DO NOT DO THIS IN A PUBLICATION***. 

It should be fairly obvious why you don't want to do this.  For instance, 
- The results (at least for the intercept) are not in agreement, given the uncertainties on the parameters
- How would you combine these results to give some meaningful answer?
- Do you actually trust the uncertainties on these parameters?  Are the parameter uncertainties actually symmetric?
- Most importantly, neither method actually accounts for both uncertainties or their correlation coefficients ($\rho_{xy}$).

Least-squares fitting reduces to some linear algebra. [See the Hogg, Bovy, & Lang paper referenced above for the derivation of these relations.]

Below we demonstrate the same results as `np.polyfit()` above.


```python
#a little bit more thorough and returns the same result as the "forward" fit
Y = obs_y.reshape(-1,1)
A = np.vstack((np.ones_like(obs_x), obs_x)).T
C = np.diag(sigma_y**2)

X = np.linalg.inv(A.transpose()@np.linalg.inv(C)@A) @ (A.transpose()@np.linalg.inv(C)@Y)

best_fit = np.poly1d(X[::-1,0])
print("The best-fit value for the slope and intercept are: {:.4f} and {:.4f}".format(X[1][0], X[0][0]))
```

OK.  The naive first attempt is not appropriate for this problem.  That least squares fitting makes certain assumptions (e.g., that the $x$ uncertainties are negligible).  Let's relax those assumptions and frame the problem more generally.

First, some statistics... (Some of you may enjoy this, others may not.  It's OK if you don't fully understand all of this today.  We will get to the `emcee` fitter very soon.)

## Minimizing $\chi^2$

A common procedure taught in undergraduate lab classes is to minimize the "chi squared" in order to identify the "best-fit" parameters $m$ and $b$: 

$$ \chi^2 = \Sigma_{i = 1}^{N} \frac{\left[y_i - f(x_i)\right]^2}{\sigma_{y_i}^2},$$

where $f(x) = mx + b$ is the model for the data. 

A least-squares fit, such as the ones used above, minimizes the value of $\chi^2$, which in turn provides a "best-fit" estimate of $m$ and $b$.

## Maximizing the Likelihood

Alternatively, we could consider a generative model, in which case the probability of any observation $y_i$ is dependent on $x_i$, our model $f(x)$, and some scatter associated with the measurement itself $\sigma_{yi}$:

$$ p\left(y_i\, | \, x_i, \sigma_{yi}, f(x_i) \right) = \frac{1}{\sqrt{2\pi\sigma_{y_i}^2}}\,\exp\left(-\frac{\left[y_i - f(x_i)\right]^2}{2\,\sigma_{y_i}^2} \right).$$

An alternative to $\chi^2$ minimization is to maximize the probability of the observed data given the model, which is to say we want to maximize the *likelihood of the observed data given the model*. We can write the likelihood function $\mathscr{L}$ as the product of the conditional probability of every observation:

$$ \mathscr{L} = \prod_{i=1}^{N} p\left(y_i\, |\, x_i, \sigma_i, f(x_i) \right)$$

This product is often difficult to evaluate, so we often take the logarithm of the likelihood, and maximize that. In this case we get:

$$ \ln \mathscr{L} = K - \sum_{i=1}^{N} \frac{\left[y_i - f(x_i) \right]^2.}{2\sigma_{yi}^2} 
   = K - \frac{1}{2}\chi^2 $$
   
where $K$ is some constant. As you can see, maximizing the likelihood (in this case) is equivalent to minimizing $\chi^2$.

Note: an important assumption in all of this is that the uncertainties, $\left(\sigma_x, \sigma_y\right)$, are Gaussian.  This is usually a reasonable assumption.

When maximizing the likelihood we can get a "pretty good" estimate of the model parameters, $m$ and $b$. This estimate is only "pretty good" because Thomas Bayes is rolling in his grave. Bayes' theorem tells us about conditional probability distributions:

$$ p\left(\theta\, |\, X, I\right) = \frac{p\left(X\, |\, \theta, I\right) p\left(\theta\, |\, I\right)}{p\left(X\, |\, I\right)} $$

In words, we want to know the (*posterior*) probability distribution for the parameters $\theta$, which $= (m,b)$ in this case, given the data $X$ and any prior knowledge $I$, $p\left(\theta\, |\, X, I\right)$.  

$p\left(X\, |\, \theta, I\right)$ is the *likelihood*, and has been defined above.  

$p\left(\theta\, |\, I\right)$ is the "*prior*" probability distribution for the parameters, that contains all of our knowledge of what those parameters should be (this can come from physical models or previous experiments).  

The denominator can be thought of as a normalization constant that most people ignore (and we will also ignore it).

Ok.  That's very nice.  

But how do we actually use all this math??!!  

And, hey, you said this was a workshop about emcee!


```python
# compute the log likelihood
def lnlike(theta, x, y, yerr):
    m, b = theta
    model = m * x + b
    ln_l = np.sum( np.log( 1./np.sqrt(2*np.pi *sigma_y**2.) * np.exp(-0.5*((y - model)/sigma_y)**2) ) )
    
    return ln_l.squeeze()

# compute the log prior
def lnprior(theta):
    m, b = theta
    if ( (0 < m < 10) and (-200 < b < 200) ):
        return 0.0
    return -np.inf

# compute the log of the likelihood multiplied by the prior
def lnprob(theta, x, y, yerr):
    lnp = lnprior(theta)
    lnl = lnlike(theta,  x, y, yerr)
    if (not np.isfinite(lnp) or not np.isfinite(lnl)):
        return -np.inf
    return lnp + lnl

```

Now, as a check, let's try to maximize this likelihood with scipy and confirm that we get the same results as the least-squares procedure described above (actually we'll minimize the negative log likelihood as this is an easier computational task).


```python
def ilnlike(theta, x, y, yerr):
    lnl = -lnlike(theta, x, y, yerr)
    return lnl.squeeze()

guess0 = np.array([2.*np.random.random() * mFreq, 2.*np.random.random() * bFreq])
scFit = minimize(ilnlike, guess0, args=(obs_x, obs_y, sigma_y))
print(scFit.message)
print("The initial guess for the slope and intercept are: {:.4f} and {:.4f}"\
      .format(guess0[0], guess0[1]))
print("The best-fit value for the slope and intercept are: {:.4f} and {:.4f}".format(scFit.x[0], scFit.x[1]))

```

## *emcee*

Now, let's use *[emcee](http://dfm.io/emcee/current/)*, a pure-Python implementation of [Goodman & Weare’s Affine Invariant Markov chain Monte Carlo (MCMC) Ensemble sampler](https://msp.org/camcos/2010/5-1/p04.xhtml) written by [Dan Foreman-Mackey](http://dfm.io/).

We've already spent enough time on background, so I will let you read about MCMC on your own (e.g., [here's the wikipedia entry](https://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo)).  You may also want to read [this journal article about emcee](https://arxiv.org/abs/1202.3665).

In short, we can use *emcee* to sample the posterior distribution that is defined by a Bayesian likelihood and priors (using MCMC).  *emcee* will use a number of "walkers" to explore this parameter space, each sampling around the maximum of the likelihood function (while accounting for the priors).  The combination of the paths that each of these walkers take, i.e. their "chains", define the posterior distribution, and provide us with probability distributions for each of the model parameters.


```python
nwalkers = 100
nfac = [1e-2, 1e-2]
ndim = len(guess0)
pos = [guess0 + nfac * np.random.randn(ndim) for i in range(nwalkers)]

sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob, args=(obs_x, obs_y, sigma_y), threads = ncores)

nsamples = 2000
foo = sampler.run_mcmc(pos, nsamples)

```

And now let's explore the output


```python
# set a "burn-in" limit
nburn = 200
m_samples = np.array(sampler.chain[:,nburn:,0]).flatten()
b_samples = np.array(sampler.chain[:,nburn:,1]).flatten()

#print the results taking 1-sigma confidence intervals
samples = np.vstack([m_samples, b_samples]).T
m_mcmc, b_mcmc = map(lambda v: (v[1], v[2]-v[1], v[1]-v[0]), zip(*np.percentile(samples, [16, 50, 84], axis=0)))
print("emcee results with 1-sigma uncertainties\n \
      m = {:.4f} +{:.4f} -{:.4f}\n \
      b = {:.4f} +{:.4f} -{:.4f}\n" \
      .format(m_mcmc[0], m_mcmc[1], m_mcmc[2], b_mcmc[0], b_mcmc[1], b_mcmc[2]))
```


```python
#some function to make plots (will re-use them later)
#to plot the chains
def plotChains(sampler, nburn, paramsNames):
    Nparams = len(paramsNames)
    fig, ax = plt.subplots(Nparams,1, figsize = (8,2*Nparams), sharex = True)
    fig.subplots_adjust(hspace = 0)
    ax[0].set_title('Chains')
    xplot = range(len(sampler.chain[0,:,0]))

    for i,p in enumerate(paramsNames):
        for w in range(sampler.chain.shape[0]):
            ax[i].plot(xplot[:nburn], sampler.chain[w,:nburn,i], color="0.5", alpha = 0.4, lw = 0.7, zorder = 1)
            ax[i].plot(xplot[nburn:], sampler.chain[w,nburn:,i], color="k", alpha = 0.4, lw = 0.7, zorder = 1)
            
            ax[i].set_ylabel(p)
            
    return ax

paramsNames = ['m','b']

axC = plotChains(sampler, nburn, paramsNames)
```


```python
#corner plot
def makeCorner(sampler, nburn, paramsNames, quantiles=[0.16, 0.5, 0.84]):
    samples = sampler.chain[:, nburn:, :].reshape((-1, len(paramsNames)))
    f = corner.corner(samples, labels = paramsNames, quantiles = quantiles)

makeCorner(sampler, nburn, paramsNames)
```


```python
#Adam's plots
def prettyPlot(x, y, xerr, yerr, rhoxy, m_samples, b_samples, m, b, nburn, ndraw = 20, extent = None):
    fig, ax = plt.subplots(1,2, figsize = (8,4))
    ax[0].set_xlabel("m")
    ax[0].set_ylabel("b")
    ax[1].set_xlabel("x")    
    ax[1].set_ylabel("y")

    ax[0].hexbin(m_samples[nburn:], b_samples[nburn:], gridsize = 250,  mincnt = 1, bins = "log", extent = extent)

    plotData(ax[1], x, y, xerr, yerr, rhoxy)
    ax[1].plot([0,300], [b + 0*m, b + 300*m], lw = 2, zorder = 2)
    if (ndraw > 0):
        for rand_draw in np.random.randint(nburn, len(m_samples), ndraw):
            m_draw = m_samples[rand_draw]
            b_draw = b_samples[rand_draw]
            ax[1].plot([-1e10,1e10], [b_draw + -1e10*m_draw, b_draw + 1e10*m_draw],
                     color = "0.5", alpha = 0.4, lw = 0.7, zorder = 1)

    return ax


axA = prettyPlot(obs_x, obs_y, sigma_x, sigma_y, rho_xy, m_samples, b_samples, m_mcmc[0], b_mcmc[0], nburn, ndraw = 20)
axA[1].set_xlim(0,300)
axA[1].set_ylim(0,700)
```

OK.  But... didn't we just get the same result as before? 

And, didn't you start this workshop talking about how we need to account for the uncertainties on $x$ as well?

## Uncertainties on both axes

To also account for the uncertainties on $x$, we need do define a new likelihood function.  We've gone through a lot of statistics already, so, for the sake of time, we'll breeze through this.  A more complete derivation is included in the Hogg, Bovy & Lang paper. 

If we assume the uncertainties are Gaussian (as we are doing), we can write the probability for obtaining a measurement, 

$$ Z_i = \left[ {\begin{array}{c}
        x_i \\
        y_i \\
        \end{array}
       }
       \right] 
$$

when the "true value" (the value you would have this data point if it had been observed with negligible noise) is:

$$ Z = \left[ {\begin{array}{c}
        x \\
        y \\
        \end{array}
       }
       \right] 
$$

as 

$$ p\left(x_i, y_i\, |\, S_i, x,y\right) = \frac{1}{2\pi \sqrt{\det(S_i)}} \exp\left(-\frac{1}{2}\left[Z_i - Z\right]^T S_i^{-1} \left[Z_i - Z\right]\right).$$

But how do we use this to fit a line?  (What should we use for Z?)  The "trick" is to project the 2D uncertainties onto the sub-space that is orthogonal to the line and to evaluate the projected displacements.  This will involve a change in variables from $(m, b)$ to $(\theta, b_\perp)$, as follows:

$$
\begin{align}
\theta &= \arctan m \\
b_\perp &= b \cos \theta \\
\hat{v} &= \left[ {\begin{array}{c}
        -\sin \theta \\
        \cos \theta \\
        \end{array}
       }
       \right] \\
\Delta_i &= \hat{v}\, Z_i - b_\perp \\
\Sigma_i^2 &= \hat{v}^T\, S_i\, \hat{v} \\
\ln \mathscr{L} &= K - \sum_{i=1}^{N} \frac{\Delta_i^2}{2\Sigma_i^2} \\
\end{align}
$$

In short, $\Delta_i$ is the orthogonal displacement of each data point $(x_i, y_i)$from the line, and $\Sigma^2$ is the orthogonal variance (by projecting the covariance matrix $S_i$, which we defined previously).

Whew!  Let's see that in action.


```python
def lnlike2(theta, x, y, xerr, yerr, rhoxy):
    th, bperp = theta
    lnl = 0.
    v = np.array([[-np.sin(th)], [np.cos(th)]])
    vT = v.transpose()
    for (xx, yy, sx, sy, rxy) in zip(x, y, xerr, yerr, rhoxy):
        S = np.array([[ sx**2, rxy*sx*sy],
                    [rxy*sx*sy, sy**2.]]) 
        Z = np.array([[xx],[yy]])
        Delta = vT @ Z - bperp
        Sigma2 = vT @ S @ v
        lnl -= Delta**2. / (2. * Sigma2)      #this ignores K (a constant), which doesn't change the likelihood maximum
        
    return lnl#.squeeze()

def lnprior2(theta):
    th, bperp = theta
    if ( (0 < np.abs(th) < 2.*np.pi) and (-200 < bperp < 200)):
        return 0.0
    return -np.inf

def lnprob2(theta, x, y, xerr, yerr, rhoxy):
    lnp = lnprior2(theta)
    lnl = lnlike2(theta,  x, y, xerr, yerr, rhoxy)
    if (not np.isfinite(lnp) or not np.isfinite(lnl)):
        return -np.inf
    return lnp + lnl

def thFromM(m):
    return np.arctan(m)
def bpFromBth(b, th):
    return b * np.cos(th)
def mFromTh(th):
    return np.tan(th)
def bFromThBp(th, bp):
    return bp/np.cos(th)
```

## Exercise:

Run this probability distribution through *emcee*, then print and plot the results.



```python
#Hint: find the related cell(s) above; copy the contents and paste it here.  Then make any necessary edits.
```


```python
# set a "burn-in" limit
nburn = 200
th_samples = np.array(sampler.chain[:,nburn:,0]).flatten()
bp_samples = np.array(sampler.chain[:,nburn:,1]).flatten()
m_samples = mFromTh(th_samples)
b_samples = bFromThBp(th_samples, bp_samples)

#print the results taking 1-sigma confidence intervals
samples = np.vstack([m_samples, b_samples]).T
m_mcmc, b_mcmc = map(lambda v: (v[1], v[2]-v[1], v[1]-v[0]), zip(*np.percentile(samples, [16, 50, 84], axis=0)))
print("emcee results with 1-sigma uncertainties\n \
      m = {:.4f} +{:.4f} -{:.4f}\n \
      b = {:.4f} +{:.4f} -{:.4f}\n" \
      .format(m_mcmc[0], m_mcmc[1], m_mcmc[2], b_mcmc[0], b_mcmc[1], b_mcmc[2]))
```


```python
paramsNames = [r'$\theta$', r'$b_\perp$']

#chain plot
axC = plotChains(sampler, nburn, paramsNames)
```


```python
#corner plot
makeCorner(sampler, nburn, paramsNames)
```


```python
#Adam's plots
axA = prettyPlot(obs_x, obs_y, sigma_x, sigma_y, rho_xy, m_samples, b_samples, m_mcmc[0], b_mcmc[0], nburn, ndraw = 20)
axA[1].set_xlim(0,300)
axA[1].set_ylim(0,700)
```

So far we haven't really showcased the power of Bayesian model fitting.  One example of when you would really want to use a Baysian tool is when you have model parameters that you want to *marginalize over*.  The following examples dive into this. 

We will first increase the scatter of our data by adding a few outliers.


## Instrinsic Scatter

Let's assume that the data is still best described by a line, and we think the best fit solution is what we found above (without the extra scatter).  


```python
obs_x = np.append(obs_x, np.array([201., 244., 287.]))
obs_y = np.append(obs_y, np.array([592., 401., 402.]))
sigma_x = np.append(sigma_x, np.array([9., 4., 7.])) 
sigma_y = np.append(sigma_y, np.array([61., 25., 15.]))
rho_xy = np.append(rho_xy, np.array([-0.84, 0.31, -0.27]))
```

## Exercise:

Use our previous likelihood function with scipy.minimize to fit these data.  Print and plot the results.



```python
#Hint: find the related cell(s) above; copy the contents and paste it here.  Then make any necessary edits.
```

The fit is different, clearly pulled by the outliers.  Both fit parameters are $>3\sigma$ away from what we found previously.   Also, notice that (at least on my notebook) the minimize function is giving an ominous warning.  (I would not want to publish the result of that fit!)

But we can deal with this using a slightly different likelihood function, and using *emcee*.  We will introduce one additional variable, $V$ : Gaussian variance orthogonal to the line.  Hogg et al. give us the likelihood:

$$
\ln \mathscr{L} = K - \sum_{i=1}^{N}\frac{1}{2}\ln\left(\Sigma_i^2 + V\right) - \sum_{i=1}^{N} \frac{\Delta_i^2}{2\left[\Sigma_i^2 + V\right]} 
$$

## Exercise:

Code this up into a new likelihood, prior and posterior probability distribution.  Run this through *emcee*.  Print and plot the results.  In your plot that includes the fit line(s), also show the "$1\sigma$" range in the scatter. 



```python
#Hint: find the right cell above; copy the contents and paste it here.  Then make any necessary edits.
```


```python
nwalkers = 100
nfac = [1e-2, 1e-2, 1.]
th = thFromM(mFreq)
bp = bpFromBth(bFreq, th)
guess0 = np.array([th, bp, 30.])
ndim = len(guess0)
pos = [guess0 + nfac * np.random.randn(ndim) for i in range(nwalkers)]

sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprobS, args=(obs_x, obs_y, sigma_x, sigma_y, rho_xy), threads = ncores)

nsamples = 2000
foo = sampler.run_mcmc(pos, nsamples)
```


```python
# set a "burn-in" limit
nburn = 200
th_samples = np.array(sampler.chain[:,nburn:,0]).flatten()
bp_samples = np.array(sampler.chain[:,nburn:,1]).flatten()
s_samples = np.array(sampler.chain[:,nburn:,2]).flatten()
m_samples = mFromTh(th_samples)
b_samples = bFromThBp(th_samples, bp_samples)

#print the results taking 1-sigma confidence intervals
samples = np.vstack([m_samples, b_samples, s_samples]).T
m_mcmc, b_mcmc, s_mcmc = map(lambda v: (v[1], v[2]-v[1], v[1]-v[0]), zip(*np.percentile(samples, [16, 50, 84], axis=0)))
print("emcee results with 1-sigma uncertainties\n \
      m = {:.4f} +{:.4f} -{:.4f}\n \
      b = {:.4f} +{:.4f} -{:.4f}\n \
      s = {:.4f} +{:.4f} -{:.4f}\n" \
      .format(m_mcmc[0], m_mcmc[1], m_mcmc[2], b_mcmc[0], b_mcmc[1], b_mcmc[2], s_mcmc[0], s_mcmc[1], s_mcmc[2]))
```


```python
paramsNames = [r'$\theta$', r'$b_\perp$', r'$\sqrt{V}$']

#chain plot
axC = plotChains(sampler, nburn, paramsNames)
```


```python
#corner plot
makeCorner(sampler, nburn, paramsNames)
```


```python
#Adam's plots
fac = 10.
axA = prettyPlot(obs_x, obs_y, sigma_x, sigma_y, rho_xy, m_samples, b_samples, m_mcmc[0], b_mcmc[0], nburn, ndraw = 0, \
                 extent = [m_mcmc[0] - fac*m_mcmc[2], m_mcmc[0] + fac*m_mcmc[1], b_mcmc[0] - fac*b_mcmc[2], b_mcmc[0] + fac*b_mcmc[1]])
axA[0].set_xlim(m_mcmc[0] - fac*m_mcmc[2], m_mcmc[0] + fac*m_mcmc[1])
axA[0].set_ylim(b_mcmc[0] - fac*b_mcmc[2], b_mcmc[0] + fac*b_mcmc[1])
axA[1].set_xlim(0,300)
axA[1].set_ylim(0,700)

smax = s_mcmc[0] + s_mcmc[1]
x0 = -b_mcmc[0]/m_mcmc[0] #x intercept
beta = np.arctan2(x0, b_mcmc[0]) #angle between line and y axis
alpha = np.pi/2. - beta #angle between perpendicular offset (s) and line
d = smax / np.cos(alpha) #distance to perpendicular offset in y direction
axA[1].plot([-1e10,1e10], [d + b_mcmc[0] + -1e10*m_mcmc[0], d + b_mcmc[0] + 1e10*m_mcmc[0]], ':k')
axA[1].plot([-1e10,1e10], [-d + b_mcmc[0] + -1e10*m_mcmc[0], -d + b_mcmc[0] + 1e10*m_mcmc[0]], ':k')
```

## Pruning Outliers

We will now append an additional outlier to our data.

Again, let's assume that the data is still best described by a line, and we think the best fit solution is what we found above (without the extra scatter or outliers). 


```python
obs_x = np.append(obs_x, 47.)
obs_y = np.append(obs_y, 583.)
sigma_x = np.append(sigma_x, 11.) 
sigma_y = np.append(sigma_y, 38.)
rho_xy = np.append(rho_xy, 0.64)
```

## Exercise:

Use our previous likelihood function with scipy.minimize to fit these data.  Print and plot the results.



```python
#Hint: find the related cell(s) above; copy the contents and paste it here.  Then make any necessary edits.
```

Clearly this result is different from what we found before, and again we see the ominous warning from scipy.

Given just the data, and thus without any direct knowledge which observations are "outliers", how could we verify this fit and deal with what we may think "by eye" are outliers?

An ad hoc solution would be to try some type of sigma clipping algorithm.  This is very popular, but also rather subjective.  We don't want to introduce unnecesary biases into our fit.  So we want to ***avoid sigma clipping***.  

Instead we will alter the likelihood function to describe a gaussian mixture model where some of the data that we have observed come from a distribution described by a line and the remaining observations come from a gaussian model. Ultimately, we don't care about which points are outliers, we are primarily interested in determining the values of $m$ and $b$. This is where we will truly take advantage of the Bayesian framework that we have been developing as we are going to *marginalize* over the outlier parameters to determine confidence regions on the values of $m$ and $b$.

With math: Recall our likelihood function that describes the fit line with uncertainties on both axes:

$$
\ln \mathscr{L} = K - \sum_{i=1}^{N} \frac{\Delta_i^2}{2\Sigma_i^2} 
$$

This is the "good" model that we care about.  Let's assume that some fraction of our data, $p_g$, comes from this model.  For this step, we need the constant, $K$, to be defined explicitly. So let's be a bit more thorough:

$$
\mathscr{L}_g = \sum_{i=1}^{N} \left( \frac{p_g}{\sqrt{2\pi\Sigma_i^2}} exp\left(  - \frac{\Delta_i^2}{2\Sigma_i^2}\right) \right)
$$ 

We now assume that the remaining fraction, $1-p_g$, of our data comes from a "bad" model that we *don't* care about, but is described by some variance, $V_b$, and mean value, $\Delta_b$, relative to the line:

$$
\mathscr{L}_b = \sum_{i=1}^{N} \left( \frac{1 - p_g}{\sqrt{2\pi\left(\Sigma_i^2 + V_b\right)}} exp\left(  - \frac{\left(\Delta_i - \Delta_b\right)^2}{2\left(\Sigma_i^2 + V_b\right)}\right) \right)
$$

The total log likelihood is the sum of these two log likelihoods written above"


$$
\ln \mathscr{L} = \ln \left(\mathscr{L}_g + \mathscr{L}_b \right)
$$

In principle we could also use this model to evaluate if any given data point is drawn from the outlier distribution ("bad") or the line ("good").  There's an interesting discussion on this on [Dan Foreman-Mackey's blog](http://dfm.io/posts/mixture-models/).


```python
#writing things slightly differently here so that I can return both "blobs" from emcee
def lnlikeOg(theta, x, y, xerr, yerr, rhoxy):
    th, bperp, _, _, _ = theta
    lnlg = np.array([])
    v = np.array([[-np.sin(th)], [np.cos(th)]])
    vT = v.transpose()
    for (xx, yy, sx, sy, rxy) in zip(x, y, xerr, yerr, rhoxy):
        S = np.array([[ sx**2, rxy*sx*sy],
                    [rxy*sx*sy, sy**2.]]) 
        Z = np.array([[xx],[yy]])
        Delta = vT @ Z - bperp
        Sigma2 = vT @ S @ v
        
        lnlg = np.append(lnlg, -np.log(np.sqrt(2*np.pi*Sigma2)) - Delta**2/(2*Sigma2))

    return lnlg

def lnlikeOb(theta, x, y, xerr, yerr, rhoxy):
    th, bperp, p_g, sigma_b, delta_b = theta
    sigma_b = np.clip(sigma_b, 0, 2000)         # prevent errors in likelihood eval
    lnlb = np.array([])
    v = np.array([[-np.sin(th)], [np.cos(th)]])
    vT = v.transpose()
    for (xx, yy, sx, sy, rxy) in zip(x, y, xerr, yerr, rhoxy):
        S = np.array([[ sx**2, rxy*sx*sy],
                    [rxy*sx*sy, sy**2.]]) 
        Z = np.array([[xx],[yy]])
        Delta = vT @ Z - bperp
        Sigma2 = vT @ S @ v
        
        lnlb = np.append(lnlb, -np.log(np.sqrt(2*np.pi*(Sigma2 + sigma_b**2.))) - (Delta - delta_b)**2/(2*(Sigma2 + sigma_b**2.)))        

    return lnlb

#priors
def lnpriorO(theta):
    th, bperp, p_g, sigma_b, delta_b = theta
    if ( (0 < np.abs(th) < 2.*np.pi) and (-200 < bperp < 200) and (0 <= sigma_b < 2000)
        and (0 < p_g < 1) and (-2000 < delta_b < 2000)):
        return 0.0
    return -np.inf


def lnprobO(theta, x, y, xerr, yerr, rhoxy):
    th, bperp, p_g, sigma_b, delta_b = theta

    lnp = lnpriorO(theta)  
    p_g = np.clip(p_g, 1e-20, 1 - 1e-20)         # prevent errors in likelihood eval

    lnlg = lnlikeOg(theta,  x, y, xerr, yerr, rhoxy) 
    arg1 = lnlg + np.log(p_g)
                           
    lnlb = lnlikeOb(theta,  x, y, xerr, yerr, rhoxy)
    arg2 = lnlb + np.log(1. - p_g)
                              
    lnl = np.sum(np.logaddexp(arg1, arg2))
                      
    if (not np.isfinite(lnp) or not np.isfinite(lnl)):
        return -np.inf, None
    
    # We're using emcee's "blobs" feature in order to keep track of the
    # "good" and "bad" likelihoods for reasons that will become
    # clear soon.
    
    return lnp + lnl, (arg1, arg2)
```

## Exercise:

Run this probability distribution through *emcee*, then print and plot the results. (Hint, you may need to increase the number of samples and the burn-in limit because we are now fitting to more parameters.)


```python
#Hint: find the related cell(s) above; copy the contents and paste it here.  Then make any necessary edits.
```

# Now let's work with real data

I pulled Table 1 from [this paper on the "M-$\sigma$" relation](http://adsabs.harvard.edu/abs/2009ApJ...698..198G)

I will use *pandas* to read in and reformat the data, but feel free to do this instead with your own favorite method.


```python
import pandas as pd
from IPython.display import display
```


```python
#data from : http://adsabs.harvard.edu/abs/2009ApJ...698..198G
#They find log(MBH/M) = α + β log(σ/200 km s−1) with (α,β,0) = (8.12 ± 0.08, 4.24 ± 0.41, 0.44 ± 0.06) 
#But they include some with upper limits that I don't have here
mSigma = pd.read_csv('apj303542t1_ascii.txt', quotechar='"', sep='\t')
#print(mSigma.columns.values)
display(mSigma)

def fixSci(series):
    x = series.str.split(' ', expand = True)
    x2 = x[:][2].str.split('^', expand = True)
    return(x[:][0].astype("float") * 10.**x2[:][1].astype("float"))

def fixSig(series):
    x = series.str.split(' ', expand = True)
    return(x[:][0].astype("float"))
def fixSige(series):
    x = series.str.split(' ', expand = True)
    x2 = x[:][2].str.split('^', expand = True)
    return(x2[:][0].astype("float"))


obs_x = np.array(fixSig(mSigma['sigma_e (km s^-1)']))
obs_lx = np.log10(obs_x)
sigma_lx = 1./(obs_x * np.log(10.)) * np.array(fixSige(mSigma['sigma_e (km s^-1)']))

obs_y = np.array(fixSci(mSigma['M_BH (M_sun)']))
obs_ly = np.log10(obs_y)

#this is probably not what we would want for publication, but we can first try to take a mean uncertainty in x
y_low = np.array(fixSci(mSigma['M_low (M_sun)']))
y_high = np.array(fixSci(mSigma['M_high (M_sun)']))
ye_low = obs_y - y_low
ye_high = y_high - obs_y
ye = np.array([ (y1 + y2)/2. for (y1, y2) in zip(ye_low, ye_high)])
sigma_ly = 1./(obs_y * np.log(10.)) * ye

rho_lxy = np.full_like(obs_x, 0.0)


```

## Exercise:

Plot these data ("$x$" is $\log_{10}\sigma$,  and "$y$" is $\log_{10} M$).  Then run the data through *emcee*, using the likelihood that accounts for intrinsic scatter.  Print and plot the results.


```python
#Hint: find the related cell(s) above; copy the contents and paste it here.  Then make any necessary edits.
```

f4e8aefe7a2a94cf9c0dcf87425f8cb375334c1c
ipynb
Jupyter Notebook
emceeDemoCRD2018_EXERCISES.ipynb
adamamiller/CRD2018emcee
3fc256f46fe89dff1a721cb4477191828af0526c
MIT
2021-03-11T18:41:26.000Z
2021-03-11T18:41:26.000Z
emceeDemoCRD2018_EXERCISES.ipynb
steveschulze/CRD2018emcee
3fc256f46fe89dff1a721cb4477191828af0526c
MIT
emceeDemoCRD2018_EXERCISES.ipynb
steveschulze/CRD2018emcee
3fc256f46fe89dff1a721cb4477191828af0526c
MIT
2019-02-27T19:33:09.000Z
2021-03-11T18:41:24.000Z
Qwen/Qwen-72B
1. YES
2. YES