<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 MCMC Demonstration
Markov Chain Monte Carlo is a useful technique for fitting models to data and obtaining estimates for the uncertainties of the model parameters.
There are a slew of python modules and interfaces to do MCMC including
Step1:  Emcee has multithreadding support. Set this to the number of cores you would like to use. In this demo we will use the python multiprocessing module support built in to emcee. Emcee can also use MPI if you're working on a cluster and want to distribute the job across nodes. See the documentation for that.
Step2:  Fitting a Line
Generate and Plot Some Random Data
Give it both x and y errors.
Step3:  Least-Squares Fit
(ignoring the x-errors)
Step4:  Maximum likelihood
So, we need to define a likelihood function.
Step5:  What about the Errors?
This is where MCMC comes in. But we need to add some priors for the parameters and use those priors
Step6:  That took about 10 seconds on my desktop (3.4 GHz Core i7).
What is the acceptance rate? Lore has it that this should be between $0.3-0.5$.
Step7:  This acceptance rate is okay. If it is too low, the emcee documentation suggests increasing the number of walkers until the acceptance fraction is at the desired level.
Let's visualize the chains.
Step8:  It looks like the walkers have "burned in" by 50 steps, so keep only those samples after 50 steps.
Step9:  What does this look like? Let's visualize with the traditional corner plot.
I will give it the actual line parameters with the "truths" parameter, so we can see how our results compare to the actual values.
Step10:  Now let's plot a bunch of sample fits from the MCMC chain, on top of our data and other models.
Step12:  Astrophysical Example
Step13:  MCMC
This will show how you might use informative priors.
Let's make sure it knows that the dust needs to be warmer than the CMB and that the amplitude needs to be positive.
Also, "normal" galaxies have dust temperatures of ~25 K, with a dispersion of a few 2K. Let's set the prior on temperature to be a Gaussian centered at 25 K with a sigma of 2.5K.
Step14:  Because of the larger parameter space and more complex model, this will take longer to run.
Step15:  Again, look at the distribution of parameter estimates. But here, show the estimated parameters from the maximum likelihood model as the "true" values.
Step16:  The offsets between the MCMC median values and the maximum likelihood are at least partially a consequence of our chosen prior on the temperature.
<ASSISTANT_TASK:>
Python Code:
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import emcee
import corner
Explanation: MCMC Demonstration
Markov Chain Monte Carlo is a useful technique for fitting models to data and obtaining estimates for the uncertainties of the model parameters.
There are a slew of python modules and interfaces to do MCMC including:
emcee
PyMC
pymultinest
emcee is fairly straightforward to use, so this demonstration is written to use that. pymultinest is worth investigating if you have large numbers of parameters (say, > 30) and/or multi-nodal solution spaces.
Required Packages
For this demo you should have the following packages installed (in addition to standard ones like numpy, scipy, matplotlib, and astropy):
emcee
corner
Optionally, install the dust_emissivity package for the last section, on fitting a blackbody.
Demo Overview
This demo will proceed very simply, and will follow the emcee tutorial for fititng a line. At the end is a short, astrophysical example, which includes non-flat priors.
Preliminaries
End of explanation
nthreads = 3
Explanation: Emcee has multithreadding support. Set this to the number of cores you would like to use. In this demo we will use the python multiprocessing module support built in to emcee. Emcee can also use MPI if you're working on a cluster and want to distribute the job across nodes. See the documentation for that.
End of explanation
# define our true relation
m_true = 1.7
b_true = 2.7
f_true = 0.3
# generate some data
N = 30
x = np.sort(10*np.random.rand(N))
yerr = 0.2+0.6*np.random.rand(N)
y = m_true*x+b_true
y += np.abs(f_true*y) * np.random.randn(N)
y += yerr * np.random.randn(N)
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.errorbar(x, y,
            yerr=yerr,
            ls='',
            marker='.',
            color='gray',
            label='Data')
ax.plot(x, m_true*x + b_true,
        color='black',
        ls='-',
        label='True Relation')
ax.set_ylabel('y', fontsize='x-large')
ax.set_xlabel('x', fontsize='x-large')
ax.minorticks_on()
ax.legend(loc='best')
Explanation: Fitting a Line
Generate and Plot Some Random Data
Give it both x and y errors.
End of explanation
A = np.vstack((np.ones_like(x), x)).T
C = np.diag(yerr * yerr)
cov = np.linalg.inv(np.dot(A.T, np.linalg.solve(C, A)))
b_ls, m_ls = np.dot(cov, np.dot(A.T, np.linalg.solve(C, y)))
print('Least squares fitting result:')
print('slope: {0:1.2f}'.format(m_ls))
print('y-intercept: {0:1.2f}'.format(b_ls))
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.errorbar(x, y,
            yerr=yerr,
            ls='',
            marker='.',
            color='gray',
            label='Data')
ax.plot(x, m_true*x + b_true,
        color='black',
        ls='-',
        label='True Relation')
ax.plot(x, m_ls * x + b_ls,
        color='red',
        ls=':',
        label='Least Squares')
ax.set_ylabel('y', fontsize='x-large')
ax.set_xlabel('x', fontsize='x-large')
ax.minorticks_on()
ax.legend(loc='best')
Explanation: Least-Squares Fit
(ignoring the x-errors)
End of explanation
import scipy.optimize as op
def lnlike(theta, x, y, yerr):
    b, m, lnf = theta
    model = m * x + b
    inv_sigma2 = 1.0/(yerr**2 + model**2*np.exp(2*lnf))
    return -0.5*(np.sum((y-model)**2*inv_sigma2 - np.log(inv_sigma2)))
# let's make some initial guesses for our parameters
# remember this is now theta and b_perp
p2 = [b_true, m_true, f_true]
nll = lambda *args: -lnlike(*args)
result = op.minimize(nll,
                     p2,
                     args=(x, y, yerr))
if not(result['success']):
    print("Max likelihood failed.")
    print(result['message'])
ml_b, ml_m, ml_f = result['x']
print("Maximum likelihood result:")
print("slope: {0:1.2f}".format(ml_m))
print("y-intercept: {0:1.2f}".format(ml_b))
print("ln(f): {0:1.2f}".format(ml_f))
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.errorbar(x, y,
            yerr=yerr,
            ls='',
            marker='.',
            color='gray',
            label='Data')
ax.plot(x, m_true*x + b_true,
        color='black',
        ls='-',
        label='True Relation')
ax.plot(x, m_ls * x + b_ls,
        color='red',
        ls=':',
        label='Least Squares')
ax.plot(x, ml_m * x + ml_b,
        color='blue',
        ls='--',
        label='Max likelihood')
ax.set_ylabel('y', fontsize='x-large')
ax.set_xlabel('x', fontsize='x-large')
ax.minorticks_on()
ax.legend(loc='best')
Explanation: Maximum likelihood
So, we need to define a likelihood function.
End of explanation
def lnprior(theta):
    b, m, lnf = theta
    if lnf >= 0.0:
        return -np.inf
    return 0.0
def lnprob(theta, x, y, yerr):
    lp = lnprior(theta)
    if not np.isfinite(lp):
        return -np.inf
    return lp + lnlike(theta, x, y, yerr)
# now let's set up the MCMC chains
ndim = 3
nwalkers = 500
steps = 500
# initialize the walkers to the vicinity of the parameters derived from
# ML
pos = [result["x"] + 1e-3*np.random.randn(ndim) for i in range(nwalkers)]
# initialze the sampler
sampler = emcee.EnsembleSampler(nwalkers,
                                ndim,
                                lnprob,
                                args=(x,
                                      y,
                                      yerr),
                                threads=nthreads)
# go! go! go!
# run the sampler for 500 steps
sampler.run_mcmc(pos, steps)
samples = sampler.chain
Explanation: What about the Errors?
This is where MCMC comes in. But we need to add some priors for the parameters and use those priors
End of explanation
print("Mean acceptance rate is: {0:1.2f}".format(np.mean(sampler.acceptance_fraction)))
Explanation: That took about 10 seconds on my desktop (3.4 GHz Core i7).
What is the acceptance rate? Lore has it that this should be between $0.3-0.5$.
End of explanation
fig = plt.figure()
dim_name = [r'$b$', r'$m$', r'$\ln f$']
for dim in range(ndim):
    ax = fig.add_subplot(ndim, 1, dim+1)
    for i in range(nwalkers):
        ax.plot(np.arange(steps),
                samples[i, :, dim],
                ls='-',
                color='black',
                alpha=10./nwalkers)
    ax.set_ylabel(dim_name[dim], fontsize='large')
ax.set_xlabel('step', fontsize='large')
Explanation: This acceptance rate is okay. If it is too low, the emcee documentation suggests increasing the number of walkers until the acceptance fraction is at the desired level.
Let's visualize the chains.
End of explanation
samples = sampler.chain[:, 50:, :].reshape((-1, ndim))
Explanation: It looks like the walkers have "burned in" by 50 steps, so keep only those samples after 50 steps.
End of explanation
fig = corner.corner(samples,
                    labels=[r"$b$", r"$m$", r"$\ln\,f$"],
                    quantiles=[0.16, 0.5, 0.84],
                    truths=[b_true,
                            m_true,
                            np.log(f_true)],
                    show_titles=True)
Explanation: What does this look like? Let's visualize with the traditional corner plot.
I will give it the actual line parameters with the "truths" parameter, so we can see how our results compare to the actual values.
End of explanation
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.errorbar(x, y,
            yerr=yerr,
            ls='',
            marker='.',
            color='gray',
            label='Data')
ax.plot(x, m_true*x + b_true,
        color='black',
        ls='-',
        label='True Relation')
ax.plot(x, m_ls * x + b_ls,
        color='red',
        ls=':',
        label='Least Squares')
ax.plot(x, ml_m * x + ml_b,
        color='blue',
        ls='--',
        label='Max likelihood')
for b, m, lnf in samples[np.random.randint(len(samples), size=100)]:
    ax.plot(x,
            m * x + b,
            color='green',
            alpha=0.1)
ax.set_ylabel('y', fontsize='x-large')
ax.set_xlabel('x', fontsize='x-large')
ax.minorticks_on()
ax.legend(loc='best')
samples[:, 2] = np.exp(samples[:, 2])
b_mcmc, m_mcmc, f_mcmc = map(lambda v: (v[1], v[2]-v[1], v[1]-v[0]),
                             zip(*np.percentile(samples, [16, 50, 84],
                                                axis=0)))
print("MCMC Parameter estimates:")
print("slope: {0:1.2f} (+{1:1.2f}, -{2:1.2f})".format(m_mcmc[0],
                                                        m_mcmc[1],
                                                        m_mcmc[2]))
print("y-intercept: {0:1.2f} (+{1:1.2f}, -{2:1.2f})".format(b_mcmc[0],
                                                        b_mcmc[1],
                                                        b_mcmc[2]))
print("\nTrue values:")
print("slope: {0:1.2f}".format(m_true))
print("y-intercept: {0:1.2f}".format(b_true))
Explanation: Now let's plot a bunch of sample fits from the MCMC chain, on top of our data and other models.
End of explanation
from dust_emissivity.blackbody import modified_blackbody
import astropy.units as u
def fit_bb(x, *p):
    
    simpler wrapper function to get the units right
    I don't care about the absolute amplitude, so the 1e-9 factor is just
    for numerical happiness.
    
    return 1.e-9* p[1] * modified_blackbody((x*u.micron).to(u.Hz, equivalencies=u.spectral()),
                              p[0] * u.K,
                              beta=p[2],
                              kappa0=0.48*u.m**2/u.kg,
                              nu0=(250*u.micron).to('Hz', u.spectral())).to('Jy').value
FIRm = np.array([(70., 50., 2.6),
                 (100., 55., 2.3),
                 (160., 34., 1.6),
                 (250., 12., 0.8),
                 (350., 4.6, 0.3),
                 (500., 1.3, 0.1)],
                dtype=[('wave', float),
                       ('flux', float),
                       ('dflux', float)])
plotrange = np.arange(FIRm['wave'][0],
                      FIRm['wave'][-1],
                      1)
def lnlike(theta, x, y, yerr):
    T, amp, beta, lnf = theta
    model = fit_bb(x, T, amp, beta)
    inv_sigma2 = 1.0 / (yerr**2 + model**2*np.exp(2*lnf))
    return -0.5 * np.sum((y-model)**2*inv_sigma2 - np.log(inv_sigma2))
# initial guesses. 25K, arbitrary 
p0 = [25, 1, 1.8, -1]
nll = lambda *args: -lnlike(*args)
maxlike = op.minimize(nll, p0,
                      args=(FIRm['wave'],
                            FIRm['flux'],
                            FIRm['dflux']),method='Nelder-Mead')
Tfit, Ampfit, betafit, lnffit = maxlike["x"]
print("Max likelihood:")
print("T: {0:1.1f} K".format(Tfit))
print("amp: {0:1.1f}".format(Ampfit))
print("beta: {0:1.2f}".format(betafit))
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.errorbar(FIRm['wave'],
            FIRm['flux'],
            yerr=FIRm['dflux'],
            ls='',
            marker='.',
            color='black',
            label='Herschel PACS+SPIRE')
ax.plot(plotrange,
        fit_bb(plotrange,
               Tfit,
               Ampfit,
               betafit),
        color='red',
        label='Max likelihood')
ax.set_ylabel(r'F$_{\nu}$')
ax.set_xlabel('$\lambda$ ($\mu m$)')
ax.set_xlim([60, 600])
ax.set_yscale('log')
ax.set_xscale('log')
ax.legend(loc='best')
Explanation: Astrophysical Example: FIR SED
Let's say we have Herschel PACS+SPIRE photometry and we want to get the dust temperature...
End of explanation
def lnprior(theta):
    T, amp, lnf, beta = theta
    if T >= 2.73 and amp > 0.:
        return -1 * (T - 25)**2 / (2 * 2.5**2)
    return -np.inf 
def lnprob(theta, x, y, yerr):
    lp = lnprior(theta)
    if not(np.isfinite(lp)):
        return -np.inf
    return lp + lnlike(theta, x, y, yerr)
ndim, nwalkers = 4, 300
pos = [maxlike["x"] + 1e-4 * np.random.randn(ndim) for i in range(nwalkers)]
Explanation: MCMC
This will show how you might use informative priors.
Let's make sure it knows that the dust needs to be warmer than the CMB and that the amplitude needs to be positive.
Also, "normal" galaxies have dust temperatures of ~25 K, with a dispersion of a few 2K. Let's set the prior on temperature to be a Gaussian centered at 25 K with a sigma of 2.5K.
End of explanation
sampler = emcee.EnsembleSampler(nwalkers, ndim, lnprob,
                                args=(FIRm['wave'],
                                      FIRm['flux'],
                                      FIRm['dflux']),
                                threads=nthreads)
sampler.run_mcmc(pos, 1000)
samples = sampler.chain[:, 100:, :].reshape((-1, ndim))
Explanation: Because of the larger parameter space and more complex model, this will take longer to run.
End of explanation
# show best-fit values as the "truth" values
fig = corner.corner(samples,
                    labels=["T", "Amp", r"$\beta$", r"$\ln\,f$"],
                    quantiles=[0.16, 0.5, 0.84],
                    show_titles=True,
                    truths=[Tfit, Ampfit, betafit, lnffit])
Explanation: Again, look at the distribution of parameter estimates. But here, show the estimated parameters from the maximum likelihood model as the "true" values.
End of explanation
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.errorbar(FIRm['wave'],
            FIRm['flux'],
            yerr=FIRm['dflux'],
            ls='',
            marker='.',
            color='black',
            label='Herschel PACS+SPIRE')
ax.plot(plotrange,
        fit_bb(plotrange,
               Tfit,
               Ampfit,
               betafit),
        color='red',
        label='Max likelihood')
for T, A, b, lnf in samples[np.random.randint(len(samples), size=100)]:
    ax.plot(plotrange,
            fit_bb(plotrange, T, A, b),
            color='green',
            alpha=0.05)
ax.set_ylabel(r'F$_{\nu}$')
ax.set_xlabel('$\lambda$ ($\mu m$)')
ax.set_xlim([60, 600])
ax.set_yscale('log')
ax.set_xscale('log')
ax.legend(loc='best')
samples[:, 3] = np.exp(samples[:, 3])
T_mcmc, A_mcmc, beta_mcmc, f_mcmc = map(lambda v: (v[1], v[2]-v[1], v[1]-v[0]),
                                       zip(*np.percentile(samples, [16, 50, 84],
                                                          axis=0)))
print("MCMC Parameter estimates:")
print("T: {0:1.2f} (+{1:1.2f}, -{2:1.2f}) K".format(T_mcmc[0],
                                                  T_mcmc[1],
                                                  T_mcmc[2]))
print("beta: {0:1.2f} (+{1:1.2f}, -{2:1.2f})".format(beta_mcmc[0],
                                                     beta_mcmc[1],
                                                     beta_mcmc[2]))
Explanation: The offsets between the MCMC median values and the maximum likelihood are at least partially a consequence of our chosen prior on the temperature.
End of explanation
<END_TASK>