1. Simple Monte Carlo Estimation Example

The probability density function (PDF) for the decay angle θ in muon decay is given by:

p(x) = 1/2 (1 + αx), where −1 ≤ x = cosθ ≤ 1, α ≈ −1/3

(a) Transformation Method:

Use the transformation method to derive an expression that allows you to generate variables x drawn from the PDF given in equation 1, using uniform random variates u drawn from U(0,1). [5 points]

(b) Simulation:

Generate 1,000 random numbers x following the PDF for α = −0.33, and plot their histogram. [5 points]

(c) Parameter Estimation:

From your simulated dataset, estimate the value of α by using:

• (i) Maximum likelihood [4 points]
• (ii) The method of moments [4 points]

2. Model Selection

The file polydat.txt (available at Canvas Link) contains 21 pairs of data points (x, y), which are most easily read using:

x, y = numpy.loadtxt(`polydat.txt`)

These data were generated using a polynomial function of order N < 10, evaluated at regularly spaced values of x, with Gaussian measurement errors μ = 0, σ = 0.1 added to each y-value.

Write a script to fit a polynomial of specified order to these data. Fit polynomials of order 0 through 9, and tabulate the following for each fit:

• • χ²

• χ² confidence level
• Akaike Information Criterion (AIC)
• Bayesian Information Criterion (BIC)

Use these values to determine the most likely polynomial order used to generate the data. Additionally, plot your polynomial fit on top of the data, with a y-error of 0.1 for each data point. [12 points]

3. MCMC Astronomy Example: Fitting the Galaxy Stellar Mass Function

In this exercise, you will:

• Write a Monte Carlo program to generate a distribution of galaxy masses following a two-parameter distribution.
• Use Markov Chain Monte Carlo (MCMC) to recover the parameters of this distribution and explore the degeneracies between them.

The stellar masses M of galaxies in the nearby Universe are described by the Schechter function:

p(log10 M | M*, α) ∝ e−M/M* (M/M*)α+1

Where α is the power-law slope at low masses, and M* represents the location of the exponential cutoff at high masses. The Galaxy and Mass Assembly Survey yields the following parameter values:

log10 M* = 10.66
α = −1.1

(a) Sampling the Distribution:

Write a program to generate the stellar masses of 1,000 galaxies in the mass range 8 < log₁₀ M < 11.5 using the distribution defined by equation (2) with the provided parameters. Use the standard U(0,1) random number generator to efficiently generate the masses. Plot a histogram of the generated masses, using 50 equally spaced bins in log mass, covering the range 8 < log₁₀ M < 11.5, with a logarithmic scale for the frequency axis. [10 points]

(b) Parameter Recovery via MCMC:

Use MCMC to recover the parameters log₁₀ M_* and α, and their 68% confidence intervals from your simulated data. Use maximum likelihood rather than a least squares fit to the histogram, and start with initial guesses of log₁₀ M_* = 11 and α = −1. Show histograms of the parameter values obtained from the MCMC chain, along with a scatter plot of log₁₀ M_* against α to illustrate any degeneracy between the parameters. It is strongly recommended that you use the example code supplied with the emcee library (available at emcee Documentation) and modify it to suit the problem, rather than writing the MCMC code from scratch. [10 points]