import numpy as np 
import matplotlib.pyplot as plt
import pymc3 as pm
import scipy.stats as stats
import matplotlib.gridspec as gridspec

fig_width_pt = 246.0                    # Get this from LaTeX using \showthe\columnwidth
inches_per_pt = 1.0/72.27               # Convert pt to inches
golden_mean = (np.sqrt(5)-1.0)/2.0         # Aesthetic ratio
fig_width = fig_width_pt*inches_per_pt  # width in inches
fig_height =fig_width*golden_mean       # height in inches
fig_size = [fig_width,fig_height]

plt.rcParams['text.latex.preamble']=[r"\usepackage{newcent}"]
plt.style.use("ggplot")
#plt.rcParams['text.usetex'] = True
plt.rcParams['font.family'] = 'serif'
plt.rcParams['figure.dpi'] = 300
plt.rcParams['figure.figsize'] = fig_size
plt.rc('xtick', labelsize='x-small')
plt.rc('ytick', labelsize='x-small')

square = [3.4039020340390205,3.4039020340390205]

Data

### As calculated in plasmoid_subband_analysis.ipynb

xdata = np.array([10.96845768, 10.97741304, 11.02116692, 11.02936182], dtype=np.float64)

y1 = np.array([2.24113874, 2.22965796, 2.1941242, 2.17601784], dtype=np.float64)
sig1 = np.array([0.01697274, 0.0176918,  0.02095898, 0.01944234], dtype=np.float64)

y2 = np.array([2.37165615, 2.36751518, 2.27429411, 2.33748569], dtype=np.float64)
sig2 = np.array([0.01709333, 0.01716664, 0.02283543, 0.01862791], dtype=np.float64)

ydata = [y1, y2]
sigma = [sig1, sig2]

Bayesian Linear Regression

! pip install emcee
 Requirement already satisfied: emcee in /home/daniel/.virtualenvs/gaston/andrew/lib/python3.6/site-packages (2.2.1)
Requirement already satisfied: numpy in /home/daniel/.virtualenvs/gaston/andrew/lib/python3.6/site-packages (from emcee) (1.16.0)
You are using pip version 18.1, however version 19.0.1 is available.
You should consider upgrading via the 'pip install --upgrade pip' command.

from __future__ import print_function, division

import os
import sys
import numpy as np

# import emcee
import emcee

# import model and data
#from createdata import *
# define the model function
def straight_line(x, m, c):
    """
    A straight line model: y = m*x + c
    
    Args:
        x (list): a set of abscissa points at which the model is defined
        m (float): the gradient of the line
        c (float): the y-intercept of the line
    """
    
    return m*x + c
def logposterior(theta, data, sigma, x):
    """
    The natural logarithm of the joint posterior.
    
    Args:
        theta (tuple): a sample containing individual parameter values
        data (list): the set of data/observations
        sigma (float): the standard deviation of the data points
        x (list): the abscissa values at which the data/model is defined
    """
    
    lp = logprior(theta) # get the prior
    
    # if the prior is not finite return a probability of zero (log probability of -inf)
    if not np.isfinite(lp):
        return -np.inf
    
    # return the likeihood times the prior (log likelihood plus the log prior)
    return lp + loglikelihood(theta, data, sigma, x)


def loglikelihood(theta, data, sigma, x):
    """
    The natural logarithm of the joint likelihood.
    
    Args:
        theta (tuple): a sample containing individual parameter values
        data (list): the set of data/observations
        sigma (float): the standard deviation of the data points
        x (list): the abscissa values at which the data/model is defined
    
    Note:
        We do not include the normalisation constants (as discussed above).
    """
    
    # unpack the model parameters from the tuple
    m, c = theta
    
    # evaluate the model (assumes that the straight_line model is defined as above)
    md = straight_line(x, m, c)
    
    # return the log likelihood
    return -0.5*np.sum(((md - data)/sigma)**2)


def logprior(theta):
    """
    The natural logarithm of the prior probability.
    
    Args:
        theta (tuple): a sample containing individual parameter values
    
    Note:
        We can ignore the normalisations of the prior here.
    """
    
    lp = 0.
    
    # unpack the model parameters from the tuple
    m, c = theta
    
    # uniform prior on c
    cmin = -8. # lower range of prior
    cmax = 30.  # upper range of prior
    
    # set prior to 1 (log prior to 0) if in the range and zero (-inf) outside the range 
    lp = 0. if cmin < c < cmax else -np.inf
    
    mmin = -5     # mean of the Gaussian prior
    mmax = 5. # standard deviation of the Gaussian prior
    lp += 0. if mmin < m < mmax else -np.inf #0.5*((m - mmu)/msigma)**2
    
    return lp
Nens = 1000   # number of ensemble points


mmin = -5     # mean of the Gaussian prior
mmax = 5.
cmin = -8.  # lower range of prior
cmax = 30.   # upper range of prior

mini = np.random.uniform(mmin, mmax, Nens)
cini = np.random.uniform(cmin, cmax, Nens) # initial c points

inisamples = np.array([mini, cini]).T # initial samples

ndims = inisamples.shape[1] # number of parameters/dimensions

Nburnin = 500   # number of burn-in samples
Nsamples = 500  # number of final posterior samples

postsamples = {}
for i in [0, 1]:

    # set additional args for the posterior (the data, the noise std. dev., and the abscissa)
    argslist = (ydata[i], sigma[i], xdata)

    # set up the sampler
    sampler = emcee.EnsembleSampler(Nens, ndims, logposterior, args=argslist)

    # pass the initial samples and total number of samples required
    sampler.run_mcmc(inisamples, Nsamples+Nburnin);

    # extract the samples (removing the burn-in)
    postsamples[i] = sampler.chain[:, Nburnin:, :].reshape((-1, ndims))

#fig.savefig('emcee.png')
# plot posterior samples (if corner.py is installed)
import corner

print('Number of posterior samples is {}'.format(postsamples[0].shape[0]))

fig = corner.corner(postsamples[0], labels=[r"$m$", r"$c$"])#, truths=[m, c])
 Number of posterior samples is 500000

png
# plot posterior samples (if corner.py is installed)
import corner

print('Number of posterior samples is {}'.format(postsamples[1].shape[0]))

fig = corner.corner(postsamples[1], labels=[r"$m$", r"$c$"])#, truths=[m, c])
 Number of posterior samples is 500000

png

Optical Thickness Diagnostic

alpha = [1.03955744367, 1.08899069766]
beta = {}
#beta[0] = [0.0, 0.405] 
#beta[1] = [0.0, 0.261]
beta[0] = [0.0, 0.458]
beta[1] = [0.0, 0.369]


def diagnostic(tau, alpha, beta):
    m = alpha*(-2*tau/(np.exp(tau)-1)) + beta
    return m 
tau = np.logspace(-3, 3, 1000)
upper = {}
lower = {}
for ib in range(2):
    upper[ib] = diagnostic(tau, alpha[0], beta[ib][1])
    lower[ib] = diagnostic(tau, alpha[1], beta[ib][0])
import matplotlib.pyplot as plt
%matplotlib inline
plt.style.use("/home/daniel/thesis/thesis-style.mpl")
ssp_legend = {'family': 'Source Code Pro',
            'weight': 'normal',
            'size': 8,
            }
#import thesis
def thesisify(f, height=1):
    from matplotlib import rc, font_manager
    lato = {'family': 'Lato',
            'color':  'black',
            'weight': 'light',
            'size': 10,
            }
    ssp_ticks = {'family': 'Source Code Pro',
            'weight': 'normal',
            'size': 6,
            }
    ssp_legend = {'family': 'Source Code Pro',
            'weight': 'normal',
            'size': 8,
            }
    ticks_font = font_manager.FontProperties(**ssp_ticks)
    # make the figure look the correct size
    #f.set_figwidth(thesis.figwidth)
    #f.set_figheight(height * thesis.figheight)
    # individual axis manipulations
    for ax in f.axes:
        for label in ax.get_xticklabels():
            label.set_fontproperties(ticks_font)
        ax.set_xlabel(ax.get_xlabel(), fontdict=lato)  
        ax.xaxis.get_offset_text().set_fontproperties(ticks_font)
        for label in ax.get_yticklabels():
            label.set_fontproperties(ticks_font)
        ax.set_ylabel(ax.get_ylabel(), fontdict=lato) 
        ax.yaxis.get_offset_text().set_fontproperties(ticks_font)
        
        if len(ax.get_ygridlines()) > 0:
            ax.grid(which="both", color='#348ABD', alpha=0.4, lw=0.3,)
        
    f.tight_layout()
    return f
cis = {}
cis[0] = np.percentile(postsamples[0][:,0], [2.5, 50, 97.5])
cis[1] = np.percentile(postsamples[1][:,0], [2.5, 50, 97.5])
for ib in range(2):
    
    
    fig = plt.figure(dpi=300)#, ax = plt.subplots(ncols = 2, nrows = 1, sharey=True)
    
    gs = gridspec.GridSpec(1, 2, width_ratios=[1, 3]) 
    
    #ax = ax.flatten()
    ax1 = plt.subplot(gs[0])
    ax = [ax1, plt.subplot(gs[1], sharey = ax1)]
    
    ax[1].fill_between(tau, upper[ib], lower[ib], alpha = 0.5, color='#2F4F4F')#556B2F', alpha=0.5)
    ax[1].set_xscale('log')
    ax[1].set_xlim(1e-3, 1e3)
    #ax[1].set_title('Box '+str(ib+1))
    ax[1].set_xlabel('Optical Thickness')
    ax[0].set_ylabel('Logarithmic Spectral Gradient')
    #density = stats.gaussian_kde(trace[ib]['gradient'])
    #ax[0].plot(density)
    ax[0].hist(postsamples[ib][:,0], bins = 40, orientation='horizontal', color = '#9f1d35',
               histtype='step', alpha = 0.6, density=True)
    ax[0].invert_xaxis()
    ax[0].set_xlabel('Occurence')
    ns = np.linspace(0,ax[0].get_xlim()[0], len(tau))
    for ci in cis:
        ax[0].fill_between(ns, np.ones(len(tau))*cis[ib][0], np.ones(len(tau))*cis[ib][2], 
                       alpha = 0.35, color = '#9f1d35')
        ax[1].fill_between(tau, np.ones(len(tau))*cis[ib][0], np.ones(len(tau))*cis[ib][2],  
                       alpha = 0.35, color = '#9f1d35')
    #ax[0].set_xlim(150, 0)   
    ax[0].set_ylim(-2.5, 0.5)
    ax[1].set_ylim(-2.5, 0.5)
    plt.setp(ax[1].get_yticklabels(), visible=False)
    plt.plot(np.ones(100), np.linspace(-2.5,0.5, 100), linestyle='--', color = '#4682b4', linewidth=0.5)
    plt.tight_layout()
    #plt.savefig('/home/arodger/Pictures/Plasmoid_project/diagnostic_box'+str(ib+1)+'.pdf')
    #plt.show()
    
thesisify(fig);
png
png
!pip freeze
 backcall==0.1.0
bleach==3.1.0
corner==2.0.1
cycler==0.10.0
decorator==4.3.0
defusedxml==0.5.0
emcee==2.2.1
entrypoints==0.3
h5py==2.9.0
ipykernel==5.1.0
ipython==7.2.0
ipython-genutils==0.2.0
ipywidgets==7.4.2
jedi==0.13.2
Jinja2==2.10
joblib==0.12.5
jsonschema==2.6.0
jupyter==1.0.0
jupyter-client==5.2.4
jupyter-console==6.0.0
jupyter-core==4.4.0
kiwisolver==1.0.1
MarkupSafe==1.1.0
matplotlib==3.0.2
mistune==0.8.4
nbconvert==5.4.0
nbformat==4.4.0
notebook==5.7.4
numpy==1.16.0
pandas==0.23.4
pandocfilters==1.4.2
parso==0.3.1
patsy==0.5.1
pexpect==4.6.0
pickleshare==0.7.5
prometheus-client==0.5.0
prompt-toolkit==2.0.7
ptyprocess==0.6.0
Pygments==2.3.1
pymc3==3.6
pyparsing==2.3.1
python-dateutil==2.7.5
pytz==2018.9
pyzmq==17.1.2
qtconsole==4.4.3
scipy==1.2.0
seaborn==0.9.0
Send2Trash==1.5.0
six==1.12.0
terminado==0.8.1
testpath==0.4.2
Theano==1.0.4
tornado==5.1.1
tqdm==4.29.1
traitlets==4.3.2
wcwidth==0.1.7
webencodings==0.5.1
widgetsnbextension==3.4.2