REMINDERS

Submitting your homework

Before doing anything,

Make a copy of this notebook in your clone of the private course repository. Edit that copy, not the version in the public repository clone. Please follow the instructions in the private repo README by locating your solutions on a path like phys366_2019/HW_week<N>/<Your Name(s)>/ and note that case matters. Take a minute to look hard at the directory structure and choose your filenames carefully: if your solutions are not in the right place they will not be graded.
Remember that you will submit your solution by pull request to the private repository.

Submitting your project abstract

This week you will submit a project abstract. As explained in the "Project Milestones" instructions:

Refine your project ideas after recieving feedback on your pitches, and continuing to brainstorm on your own. Update the pitch (which we will now arbitrarily call an abstract) to reflect any changes.

This only applies to extant projects, naturally. That is, we only need abstracts for projects that will go forward, which is not necessarily every project that was pitched (if, e.g., some projects have been consolidated after discussions with classmates).

Note that, from here on out, these milestones will be organized in the repo by project rather than by person. Therefore, push your (team's) abstract to your fork of the private repo in a folder like phys366_2019/Project_milestones/<project name>/ and submit a PR to the private repository. Like the pitch, you can write your abstract in either markdown, plain text, or jupyter notebook form. Do include your names as authors of the abstract, so we know who is working on which project.

Sending us feedback

Once you've submitted your solution, don't forget to also fill out the very quick (and required!) weekly feedback form.

Week 5 Homework

Gibbs Sampling for the Cepheid Period-Luminosity Relation

This problem continues the week's tutorial! You'll probably want to re-use the code from that notebook. You could laboriously copy code cells from one notebook to the other. Alternatively, you could export the tutorial notebook as Python code and use exec(open('mytutorial.py').read()) here. In that case, your solution should include both that raw code and the tutorial notebook (with outputs), since that stuff is part of the solution.

Preliminaries

Reproduce the setup from the tutorial somehow (see above). This includes the definitions of REMOVE_THIS_LINE and REPLACE_WITH_YOUR_SOLUTION, unless you want to see some weird behavior.

We do provide some code below, which you are welcome to use, but it assumes that variables have the same names as in the tutorial, etc.



In [ ]:

    
try:
    exec(open('Solution/Cepheids.py').read())
except IOError:
    REPLACE_WITH_YOUR_SOLUTION()

The posterior

We're now going to fit the various parameters of the (simplified) Cepheid model. If you've managed to get through the tutorial without explicitly writing down that posterior, do it now!



In [ ]:

Solution by Gibbs sampling

Yes! You will now implement a Gibbs sampler to draw from the posterior, which is conveniently fully conjugate (with the right form of priors). If you've been following so far, you know that the model now has the following free parameters:

The true values of $a_i$ and $b_i$ for each galaxy
The hyperparameters $\bar{a}$, $\bar{b}$, $\sigma_a$, and $\sigma_b$

First, work out the fully conditional posterior for each parameter

$p(a_i|...) = $
$p(b_i|...) = $
$p(\bar{a}|...) = $
$p(\sigma_{a}|...) = $
$p(\bar{b}|...) = $
$p(\sigma_{b}|...) = $

We haven't said yet what the priors on the hyperparameters are. Look up what form would make the posterior conjugate, and assign a defensible prior. Is the conjugate prior compatible with something "uninformative" (e.g. a uniform prior or the Jeffreys prior)?

For $a_i$ and $b_i$, you might need to look up the conditional form of the 2D Gaussian distribution and/or some other Gaussian identities in order to simplify. The RHS is a single, standard PDF in every case.



In [ ]:

Next, write a function to sample each of the parameters. The beautiful thing about conjugate Gibbs sampling is that the filled-in bulleted list of conditional posteriors above translates directly into the code you need to write.

Some code is provided below to organize the parameters and plot results. Use it or not, at your discretion. Caveat emptor.



In [ ]:

    
# Here's a free little class to provide access to the free parameter values, and useful subsets of them
class Parameters:
    def __init__(self):
        # Unless you have a better idea, assume the parameter vector is ordered like
        # [a_0, a_1, ..., b_0, b_1, ..., abar, sigma_a, bbar, sigma_b]
        self.all = np.zeros(2*len(ngcs)+4) # you will specify starting values later
        self.ais = self.all[0:len(ngcs)] # should be a view
        self.bis = self.all[len(ngcs):(2*len(ngcs))] # should be a view
        self.names = np.concatenate([['a_'+str(int(ngc)) for ngc in ngcs], ['b_'+str(int(ngc)) for ngc in ngcs], 
                ['abar', 'sigmaa', 'bbar', 'sigmab']])
    def a(self): return self.ais # return the whole vector
    def b(self): return self.bis # ...
    def ai(self, i): return self.ais[i] # return just one
    def bi(self, i): return self.bis[i] # ...
    def abar(self): return self.all[-4] # get the a and b distribution parameters
    def sigmaa(self): return self.all[-3] # ...
    def bbar(self): return self.all[-2] # ...
    def sigmab(self): return self.all[-1] # ...
    def set_ai(self, i, v): self.ais[i] = v # set parameters
    def set_bi(self, i, v): self.bis[i] = v # ...
    def set_abar(self, v): self.all[-4] = v
    def set_sigmaa(self, v): self.all[-3] = v
    def set_bbar(self, v): self.all[-2] = v
    def set_sigmab(self, v): self.all[-1] = v

Declare a Parameters object. Initially, all parameters are given a value of zero in the code above, which is a problem for the sigmas. The code below avoids dividing by zero, but doesn't do anything about setting an intelligent starting point. Do so.



In [ ]:

    
p = Parameters()
p.set_sigmaa(1.0)
p.set_sigmab(1.0)

Define whatever variables you need to specifiy your priors, and a function/functions to update all of the values in the an object of class Parameters.



In [ ]:

    
try:
    exec(open('Solution/sampler.py').read())
except IOError:
    REMOVE_THIS_LINE()
    def sample(current_params):
        REPLACE_WITH_YOUR_SOLUTION()
    # probably some other utility functions, and definitions that specify priors

Create an array to hold the chain, and run the sampler!



In [ ]:

    
Nsteps = 10000 # change as you see fit
chain = np.zeros((Nsteps,len(p.all)))



In [ ]:

    
%%time
for i in range(chain.shape[0]):
    sample(p)
    chain[i,:] = p.all

Plot the traces



In [ ]:

    
plt.rcParams['figure.figsize'] = (12.0, 64.0)
for j in range(chain.shape[1]):
    plt.subplot(chain.shape[1], 1, j+1)
    plt.plot(chain[:,j])
    plt.ylabel(p.names[j], fontsize=20)

Remove burn-in; do any thinning you want to do; verify convergence and robustness.



In [ ]:

Finally, look at the posterior for the hyperparameters, and offer an answer to the question posed at the top of the notebook: is the P-L relation universal according to these data, or is there evidence for additional scatter?



In [ ]:

    
# Convenient for looking at the 1D and 2D posteriors of a few parameters
corner.corner(chain[:,18:22], show_titles=True, labels=p.names[18:22]);



In [ ]:

Bonus (not had enough yet?)

Here are some ideas for digging deeper with this problem:

Constrain the distributions of $a$ and $b$ using Metropolis sampling (or some other sampler of your choice) and compare performance.
Extend the model such that the joint distribution of $a$ and $b$ is a 2D Gaussian, with a $2 \times 2$ covariance matrix describing the scatter.
Constrain the model without using the WLS trick (sampling all of the Cepheid magnitudes) using the sampler of your choice. Who knows, maybe you'll find one that that can handle the challenge gracefully.



In [ ]: