New Reading

The reading assignment for this week is Sivia chapter 4.1 and the first subsection of 4.2 (up to 4.2.1).

Also, check out this post from last year on the Jeffreys' prior for a scale parameter, such as the slope of a linear fit. If you're feeling super inspired, check out  this dense yet handy Wikipedia entry.

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

CA3 data

Part 6 of CA3 instructs you to construct a data table based on FV05. In the interest of time and uniformity, you should just use this table instead:

http://www.astro.caltech.edu/%7Ejohnjohn/astrostats/data/fv05_vf05_data.txt

The middle column indicates 0 if the star doesn't have a planet, and 1 if it does.

Parts 1-8 will be due Monday. Continue on to problems 9 if you feel particularly inspired to learn about how to deal with measurement uncertainties when doing this sort of analysis. 

Part 10 instructs you to write up your results. This is left over from last year. This year I've already encouraged you to make your results clear, either in an iPython workbook or better yet a LaTeX document. You may turn in your completed activity the same way you did with CA1 and CA2.

Class Activity 3

http://www.astro.caltech.edu/~johnjohn//astrostats/CA3.pdf

Dispositive Nulls and Detection Limits

My good friend Prof. Jason Wright has been tackling the notion of a "dispositive null" over on his professional blog. Here's his first entry, and his followup. This is good stuff and a natural extension of the notion of assigning confidence to our belief in a hypothesis.

avoid using loops!

So one of the important things to know about python is that it is an "interpreted" language, not a compiled language.  Interpreted languages (like python, IDL, and Matlab) are designed to be "higher-level" languages, and thus do not have to be compiled in the same way as C or C++, for example.  This makes some things simpler, but there are some costs.

One of the costs is that loops take a long time to execute.  See the following example:

As you see, doing operations element-by-element on an array is MUCH slower than working with whole arrays at a time.  In fact, behind the scenes, numpy is using all the power of compiled code to make these array operations lightening-fast.  But you don't need to know how this works; all you need to do is get comfortable working with whole arrays, whenever possible.  A good example of this in action is the difference between the following two ways to create arrays of random numbers:

CA1 Solution Sets

Everyone did a nice job on the Class Activities this week. Thank you for all your hard effort and careful work. Coco has placed marked-up PDF files in all of your Dropbox folders. Here are three writeups that I feel constitute the "solution set."

Solution 1 (Peter)
Solution 2 (Melodie)
Solution 3 (Adam)

In the future, please strive to have your write-up be clear enough to serve as the solution set. Science is all about communication, so the clarity of your presentation is a big part of your assessment in this class.