2D posterior plot is easy ---CA3_problem4

Debugging is annoying and we frequently run into a lot of problems when plotting 2D posterior contour.  Believe me, you are not alone.  There is a simple code I used for CA3 problem4. It is not optimized but it shows some standard steps to do the 2D posterior plots which I find useful.

# step1: read your data and do some basic check (visualization is important)
# step2: prepare a 2D -parameter mesh grid, which can be done by meshgrid function in numpy
# step3: Given the 2D parameter grid, we can calculate the probability of each datum at the every grid    coordinate, i.e. , a 2D probability for each data point. It is usually convenient to do it in log scale.
# step4:  sum the log probability of all data points, then we have the joint 2D log probability for
 the whole data set.
# step5: find the maximum joint log probability and subtract it. We do this to avoid numerical outflow problem since unnormalized  joint probability is usually a very small value.
# step6: exp(log prob) and plot a contour

The example of CA3_problem4

from numpy import *
import math
import matplotlib
import matplotlib.pyplot as plt
import asciitable
import scipy

# read the data
data      = asciitable.read('fv05_vf05_data.txt')
detection = data['comp'] #detection or not 1/0
feh       = data['metallicity'] #metallicity

# check data

n         = detection.size    #the total number of data point
print feh.min(),feh.max()     # metallicity range in the data

#choose the proper range of alpha, f>0 --> alpha>0
alpha_grid= arange(0.01,0.1,0.001)
beta_grid = arange(.0,3.0,0.01)

# generate a 2D meshgrid for alpha and beta
(alpha,beta) = meshgrid(alpha_grid,beta_grid)
logpdf = zeros((beta_grid.size,alpha_grid.size))

# loop through each data point
for i in range(n):
    f = alpha*10**(beta*feh[i])
   
    f[(f >1)] = None
    #set f to none if f>1, here we don't need to check f>0
    # because we already set alpha>0
   
    # sum the log(prob) of the all data
    logpdf += log(detection[i]*f+(1-detection[i])*(1.-f))

# set the log(prob) at forbidden (alpha,beta) value to -999   
whereAreNaNs = isnan(logpdf)
logpdf[whereAreNaNs] = -999

# find the max logpdf
print logpdf.max()

# find the peak location

ind = where(logpdf == logpdf.max())

print 'the peak is at alpha=%.4f,beta=%.4f' %(alpha_grid[ind[1]],beta_grid[ind[0]]) 

logpdf = logpdf-logpdf.max()

# pdf is the normal probability function
pdf = exp(logpdf)

# 2D plot

plt.figure()
CS = plt.contour(alpha,beta,pdf)

Avoiding Loops: a homework example

In previous post, Tim has shown us that using array is much efficient than loops. Here I'm going to give a real example about how to use numpy array to solve one of our homework problem.

The example is CA4, problem3. We have a data set of the masses of known transit planets and want to find the power index of the mass distribution function.


##########################################################

from numpy import *
import asciitable
import pylab as py

# read the transit data
all_transits = asciitable.read('transits.csv',delimiter=",")
transit_mass = all_transits['MSINI']   

# find the data with 1 < m <20
ind = ((transit_mass >= 1)& (transit_mass <=20)).nonzero()[0]
mass_data = transit_mass[ind]

n =  mass_data.size  

# the total number of planets in the given mass range

# build an array for a possible alpha range
alpha_grid = arange(1.1,3.5,0.001)

# here we use numpy array rather than loop
# (1) we build 2D meshgrid arrays for alpha and mass
# NumPy provides a convenient function for generating meshgrid arrays,
# called, appropriately, meshgrid.

(alpha,mass) = meshgrid(alpha_grid,mass_data)
#print alpha
#print mass

#(2) Calculate the probability at each alpha and m
# Attention: array math operates on an element by element basis,
# i. .e, unlike matrix math.
mass_max =  mass_data.max()
mass_min = mass_data.min()
const = 1./(1-alpha)*(mass_max**(1.-alpha)-mass_min**(1.-alpha))
log_pdf = log(1./const*mass**(-1.*alpha))

# (3) we sum the log_pdf in the mass dimension
log_pdf_alpha = log_pdf.sum(axis = 0)
plot(alpha_grid,exp(log_pdf_alpha-log_pdf_alpha.max()))

##########################################################

The example code and transit planet data can be downloaded from here:

I find this introduction of numpy calculation very useful. It is totally worth your time.

http://www.physics.nmt.edu/~raymond/software/python_notes/paper003.html

arbitrary precision in Python

Some of you mentioned about the rounding error in python. There is a Python package which can give you arbitrary precision. Take a look if you are interested.

http://code.google.com/p/mpmath/

On Python Functionality

"Any self-respecting programming language has a histogram function!" -Jennifer

"Oh! Here it is!" -Jessica

IDL to Python translation guide

https://www.cfa.harvard.edu/~jbattat/computer/python/science/idl-numpy.html

Programming Investment: Making Figures

If you continue on in astronomy, you'll find (if you haven't already) that most of the actual time you spend doing research is, in fact, programming. Not just programming to actually do whatever analysis you're doing, but often a comparable amount of time putting together nifty professional-looking, presentable figures to communicate your results.

With this in mind, you should be thinking about building an accessible toolkit of useful routines that perform tasks that you expect to repeat often. A first example of something like this would be the function we had you write a few weeks ago that calculates a 1-D confidence interval, given vectors of x and p(x). You should consider writing similar routines for common plotting tasks.

For example, I've written a function called "plot_posterior" that takes as input x and p(x) (or L(x)--not necessarily normalized), just like your confidence interval function, and makes a pretty plot, as follows:

>>> x = arange(-2,6,0.001) #make x vector

>>> px = normal(x,2,1) #evaluate N(2,1) on x

>>> plot_posterior(x,px,name='x',labelpos=(0.65,0.8)) #make the plot

(Note the default thicker-than-normal axes and lines, which is the quickest way to make your plots look more presentable and should become a habit.)

The function also has other optional parameters, such as plotting different confidence intervals, marking the median value instead of the mode as the representative value, using a symmetric interval instead of the shortest interval, turn off shading, turn off labeling, change the number of decimal places to which the result and errors are quoted, etc. This way, whenever I want to make a nice plot of a posterior distribution with regions marked and labeled and all, I don't have to start from scratch.

I've also written a similar function ('plot_posterior2d') to summarize the results of a two-dimensional parameter estimation, as you are currently working on with the Fischer & Valenti data. This produces the following, which is a nice concise way to illustrate the results:

>>> plot_posterior2d(a,b,L,confs=[0.68,0.95,0.99],'a','b',decplaces1=3)

As you can see, this function uses the 1D version of 'plot_posterior' to make the marginalized pdfs (once called with the 'horizontal' option). (It also draws only selected confidence contours on the contour plot, which is something we have not yet asked you to do, but you should be thinking about.) And now, with a modest front-end time investment, I can instantly make a snappy publication-quality plot of any two-dimensional likelihood surface I wish, potentially saving lots of time down the road.

NOTE: I'm doing all my work here using python, which I heartily recommend to all of you to experiment with at some point, if solely because of this amazing collection of example plots that I visit anytime I ever want to make any kind of fancy-looking plot. The language is quite similar in many ways to both Matlab and IDL, and has the benefits that it's both much more portable and much more free than either. I've just picked it up in the last year or so, and I'd be happy to help guide anyone who's interested.