## Monday, August 3, 2015

### Astro-Stats & Python : Bayes' Theorem (for Astrophysics)

Back to the Astro-Stats!

In my Bayes' Theorem (for a Coin) post, I introduced Bayes' Theorem for probability along with an example with a Python script.  This time, I will use the same concept for astrophysical data.  Let me first set the stage for this new example and then I'll get to the code so that you can follow along your computer.

My data comprises of the metallicity of stars and the confirmation/denial of a planet detection around a star.  You can see the text file containing this data here.  If you want to follow along, you must copy this information and paste it into your own text file; or, you can download this file using the "wget" command (on that same link) in the Terminal.

 Figure 1: Snapshot of part of table in the ".txt" file.  The first column is the real name of the star.  The second column is the metallicity of the star relative to the Sun.  This metallicity is the ratio of iron abundance to hydrogen abundance [Fe/H].  Because the Sun is used as the reference, Fe/H for the Sun equals 0.
As with my last analysis with Bayes' Theorem, I must input a prior for each parameter sought, find the Likelihood and then evaluate the new parameters.  In terms of mathematics, the differences between this exercise and the Coin example are conspicuously seen within my Likelihood and Bayes' Probability equations.

The Likelihood is equal to the probability (for the parameters' being measured) of all of the stars with detected planets multiplied by the probability (for the same parameters) of all of the stars without a planet detection.  Because each detection or absence thereof are events that occur independently of each other, the probabilities must be multiplied together to obtain the total probability for the stars with detected planets, as well as those without planets.  This can be described as
$\mathcal{L}=\prod^{k} P(\alpha,\beta) \cdot \prod^{N-k}[ 1-P(\alpha, \beta) ]$
where $$\mathcal{L}$$ is the Likelihood, k is the number of stars with planets, N is the number of stars total, P is the probability of stars with planets, and $$\alpha$$ and $$\beta$$ are parameters.  The rightmost $$\prod$$ term is the total probability of stars without detected planets and the $$\prod$$ term to the left of that is the total probability of stars with detected planets.

Essentially, the Likelihood is being found by multiplying two probabilities.  This corresponds to the generic definition of Likelihood, which can be written as $$\mathcal{L} = \prod^{N}_{i}P_{i}$$  where N is the total number of probabilities.  Likelihood can be defined as the product of all the probabilities.  This was discussed briefly in my Exercise in Likelihood and $$\chi^2$$, continued... post.

For this specific example, the given probability equation is defined as
$P (\alpha, \beta) = \alpha 10^{\beta M}$
where M is the metallicity of one of the observed stars.  Be aware that many statistical applications assume the probability of all outcomes to be equal to each other.  This equation signifies that the probability of each star for having a planet companion is distinct and depends primarily on the metallicity of that star along with the parameters.

Unlike the Coin example, I must now suggest many priors for my parameters $$\alpha$$ and $$\beta$$.  Previously, I provided only 1 prior probability and then let my function use the newfound posterior probability as the subsequent iteration's prior probability.  My methodology is different now because I will use the Grid Search and contour plot strategy to determine which of my priors yielded the highest Likelihoods and thus the highest posterior probabilities.  This method was used in my Grad Level to Intro to Astro-Stats for Undergraduates post and its sequel.  Similar to the equation stated for my Coin example, the posterior Likelihood is described as the following proportion.
$\mathcal{L}_{posterior} \propto \alpha \cdot \beta \cdot \mathcal{L}_{prior}$
Each $$\mathcal{L}_{posterior}$$ will be tracked inside of a matrix, or grid, of numbers that will later be used to create a contour plot.

Now that the math and concepts are (hopefully) clear to you, we can get to the real-life application of all of this knowledge.  My next post will illustrate the results of utilizing this statistical method in my Python script.

to be continued...