HW 5 Clarifications: What does the evidence ratio mean? (and revisions to the assignment)

Questions 4 and 5 on the HW are a bit confusing. Here's how to think of it. [Skip to the bottom if you're confused by the following ramblings about the role of Monte Carlo simulations in this kind of analysis.]

An evidence ratio (or "Bayes factor") has a specific meaning: that is, it is the odds ratio in favor of one model over another. If Ev1/Ev2 = 9, then your data give 9:1 odds that Model 1 is right; in other words, you have 90% confidence (9/(9+1)) in Model 1.

Why then, you might (rightly) ask, are we asking you do to Monte Carlo simulations? And what then is the point, you might (astutely) inquire, of Section 3.2 of our paper?

Well, to be honest, the more I think about it the more unnecessary I believe Section 3.2 is. Except for two important purposes: a sanity check and a convincing demonstration. Monte Carlo simulations allow you test the frequentist performance of your Bayesian statistic; this is what we do in Section 3.2, and what we are asking you to do in question 4. In Section 3.2 it was relatively straightforward to do Monte Carlo simulations to do this, because the two models I was comparing have no free parameters. And we find that the simulation-based confidence of 94% is (almost) exactly what you get by just calculating the confidence directly from the R value. HOWEVER, it is no longer clear to me how to do these simulations, or exactly what they mean, in cases where you have free parameters. What I did in this case (Section 4.2/4.3) was to simulate data with various values of the parameters and weight the resulting R distribution by the values of the parameter that the data favored (i.e. marginalized over the parameter). This ends up with a result that does in fact confirm the confidence that the Bayes factor suggests (R = 1.16 gives 94% confidence analytically), but this is a different procedure from what we have asked you to do in Question 4, where we tell you to simulate data just using the best fit parameters.

In addition to this, Question 5 is somewhat unclear.

John's out of town right now, but in my role as temporary deputy of this class, I ask you to combine the tasks that the current Questions 4 and 5 are asking you to do by answering the following two-part question:

How likely are you to be able to confidently (>95%) detect the presence of two lighthouses, as a function of their separation? And how likely are you to confidently select a one-LH model over a two-LH model, if the data really do come from just one LH?

This will still require you to do MC simulations, but this time not to try to confirm the frequentist performance of your statistic, but as a way to help you understand the limits and capabilities of your experiment.

To turn in Monday:

• Example 1-d and 2-d posterior plots for each data set (i.e. the one generated with 1 LH and the one generated w/ 2 LHs), with the true values of alpha marked. I would like the 2-d posteriors to have 1-, 2-, and 3-sigma contours, though if you don't get around to this yet, not a huge deal.
• A plot of "probability of confidently detecting 2 LHs" vs. LH separation (with a caption explaining what the plot means that someone not taking this class could understand)
• A histogram of R values resulting from a large number of data sets simulated using one LH, with a caption giving your interpretation of what this means.

HW4: My Solutions

Hey class- Here are the results I get from the two exercises. My prior confidence in myself doing things like this correctly on the first try is generally something around 93%. So if you get something different, check yourself first, but if your result is persistently different from mine, let me know!

Question 1: What kind of coin?

Question 2: Fair flipper or cheater?

Remember to write captions (short but to the point) for each figure: be very clear in your language about what you're calculating! This is particularly important because most people think Bayesian methods are scary and magical; your job is to show that what you do is clear and logical!

HW4, Question 2 clarification

Hey folks- I apologize for the incomplete and inexact wording of question two on the latest worksheet.

Here is the story. A Type I flipper always tells the truth. A Type II flipper, with some unknown probability f (f less than 0.5), will decide before a flip that, no matter what you guess, he will make you wrong. That is, no matter what the coin actually is, if you guess 'heads' he says 'tails,' and vice versa.

Model Selection in the Literature

Tim and I just published a paper that uses the Bayesian model selection concept we are learning about this week. The paper can be found here.

Also, here's a link to Tim's newly minted home page.

Homework 3: What is due Wednesday

For Wednesday, you can skip Problem 4. I think it is too time-intensive for the lesson I wanted you to learn. The punchline is shown in the figure below, which was generated by Jon Swift (of Swift & Beaumont 2009).
The point is that the least-squares approach underestimates the power-law index for all values of N. So the next time you see a power-law fit in the literature, be suspicious!

For Wednesday, do Problem 3. In addition to what is written, estimate the power-law index for the exoplanet masses using MLE, and compare your answer to the value you find from a least-squares fit as discussed in Problem 4. For the least-squares fit, you can ignore mass bins containing zero planets.

Turn in a figure showing the pdf for alpha, and a figure illustrating the best-fitting powerlaws from both methods compared to a histogram of the planet masses. For the two figures, write captions describing what is plotted.

ETC: less than 2 hours

Reading assignment due Wed is still Sivia Chapter 4.1

Are your likelihoods slowing you down?

Sorry about the formatting here; looks like the Google docs embed thing isn't quite compatible with the blog formatting....Hope you can still read what you need here.

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.