Part 4: Random Sampling - From Theory to Computation

Learning Outcomes¶

By the end of Part 4, you will be able to:

• Implement the inverse transform method to sample from analytical distributions including exponentials and power laws

• Apply the Cumulative Distribution Function (CDF) concept to transform uniform random numbers into any desired distribution

• Sample from complex astrophysical distributions like the Kroupa IMF using piecewise techniques

• Design rejection sampling algorithms for distributions without analytical inverse transforms

• Generate realistic initial conditions for N-body simulations including mass, position, and velocity distributions

3.1 Why Random Sampling Matters¶

Priority: 🔴 Essential Theory tells us distributions exist. But for simulations, we need to generate samples. The challenge: computers only produce uniform random numbers [0,1]. How do we transform these into complex astrophysical distributions?

This is where theory meets computation. Random sampling bridges the gap between:

• Statistical understanding (Parts 1-2)

• Computational modeling (your projects)

4.2 The Cumulative Distribution Function (CDF) and Inverse Transform Sampling¶

Priority: 🔴 Essential The key to sampling from any distribution lies in understanding the cumulative distribution function (CDF) and how it enables the inverse transform method—the most fundamental technique for converting uniform random numbers into samples from any distribution.

Understanding the CDF¶

The CDF is defined as:

$\boxed{F(x) = P(X \leq x) = \int_{-\infty}^{x} f(x') dx'}$

While the probability density function (PDF) f(x) tells us the relative likelihood at each point, the CDF F(x) tells us the accumulated probability up to that point. This accumulation is what makes sampling possible.

To understand why this is so powerful, consider what the CDF actually represents. If you have a PDF f(x) that describes the distribution of stellar masses, then F(m) answers the question: “What fraction of stars have mass less than or equal to m?”

For example, in a stellar population:

• $F(0.5 M_\odot)$ = 0.3 means 30% of stars have mass $\le 0.5 M_\odot$

• $F(1.0 M_\odot)$ = 0.6 means 60% of stars have mass $\le 1.0 M_\odot$

• $F(∞) = 1.0$ means 100% of stars have some finite mass

Key properties of CDFs:

• Always non-decreasing: $F(x₁) ≤ F(x₂) if x₁ < x₂$

• Bounded: $F(-∞) = 0$ and $F(∞) = 1$

• Continuous from the right: Important for discrete distributions

• Derivative gives PDF: $dF/dx = f(x)$ for continuous distributions

The Inverse Transform Method¶

The crucial insight is that the CDF transforms any distribution—no matter how complex—into a uniform distribution on [0,1]. This provides the bridge between uniform random numbers (which computers generate) and any distribution we want.

Visual intuition: Imagine the CDF as a transformation that “stretches” the uniform distribution. Regions where the PDF is large (high probability density) correspond to steep sections of the CDF. When we invert, these steep sections get “compressed” back, allocating more samples to high-probability regions.

Example: Exponential Distribution¶

Before tackling power laws, let’s see how this works for the exponential distribution (which describes photon path lengths, radioactive decay, and time between stellar collisions):

PDF: $f(x) = \lambda e^{-\lambda x}$ for x ≥ 0

CDF: $F(x) = \int_0^x \lambda e^{-\lambda x'} dx' = 1 - e^{-\lambda x}$

Inverse: Solve $u = 1 - e^{-\lambda x}$ for x:

• $e^{-\lambda x} = 1 - u$

• $x = -\frac{1}{\lambda}\ln(1-u)$

Since u and (1-u) are both uniform on [0,1], we can simplify to: $x = -\frac{1}{\lambda}\ln(u)$

# Sample exponential distribution (e.g., photon mean free paths)
def sample_exponential(lambda_param, n_samples):
    u = np.random.uniform(0, 1, n_samples)
    return -np.log(u) / lambda_param

# Example: mean free path of 1 pc
mfp = 1.0  # pc
lambda_param = 1.0 / mfp
path_lengths = sample_exponential(lambda_param, 10000)
print(f"Mean path: {np.mean(path_lengths):.2f} pc")
print(f"Theory: {mfp:.2f} pc")

Power Law Distributions: The Foundation of Astrophysics¶

Power law distributions appear everywhere in astronomy:

• Stellar initial mass function (IMF)

• Galaxy luminosity functions

• Cosmic ray energy spectrum

• Size distribution of asteroids and dust grains

Deriving the sampling formula:

1. Normalize the PDF: $f(x) = \frac{(\alpha-1) x^{-\alpha}}{x_{\min}^{1-\alpha} - x_{\max}^{1-\alpha}}$ (for α ≠ 1)

2. Compute the CDF: $F(x) = \int_{x_{\min}}^x f(x')dx' = \frac{x^{1-\alpha} - x_{\min}^{1-\alpha}}{x_{\max}^{1-\alpha} - x_{\min}^{1-\alpha}}$

3. Invert to get sampling formula: Set F(x) = u and solve for x: $\boxed{x = \left[x_{\min}^{1-\alpha} + u(x_{\max}^{1-\alpha} - x_{\min}^{1-\alpha})\right]^{1/(1-\alpha)}}$

Special case α = 1: The integral of 1/x is ln(x), so: $\boxed{x = x_{\min} \left(\frac{x_{\max}}{x_{\min}}\right)^u}$

Implementation and Visualization¶

The figure above shows our power law sampling in action. Here’s how to implement it yourself:

import numpy as np
import matplotlib.pyplot as plt

def sample_simple_power_law(alpha, x_min, x_max, n_samples):
    """
    Sample from a power law distribution p(x) ∝ x^(-alpha).
    
    This demonstrates the inverse transform method for power laws.
    Special case: alpha = 1 requires logarithmic sampling.
    
    Parameters:
    -----------
    alpha : float
        Power law exponent (positive for decreasing function)
    x_min, x_max : float
        Range of the distribution
    n_samples : int
        Number of samples to generate
        
    Returns:
    --------
    samples : ndarray
        Random samples from the power law
    """
    u = np.random.uniform(0, 1, n_samples)
    
    if abs(alpha - 1.0) < 1e-10:
        # Special case: p(x) ∝ x^(-1)
        # CDF integral gives logarithmic form
        samples = x_min * (x_max/x_min)**u
    else:
        # General case: p(x) ∝ x^(-alpha)
        # Apply the inverse transform formula we derived
        samples = (x_min**(1-alpha) + u*(x_max**(1-alpha) - x_min**(1-alpha)))**(1/(1-alpha))
    
    return samples

# Example: Sample stellar masses with different power laws
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# Different IMF slopes (Salpeter has alpha = 2.35)
alphas = [1.3, 2.3, 3.5]
x_min, x_max = 0.1, 100  # Solar masses

for ax, alpha in zip(axes, alphas):
    # Generate samples
    samples = sample_simple_power_law(alpha, x_min, x_max, 10000)
    
    # Plot histogram
    bins = np.logspace(np.log10(x_min), np.log10(x_max), 50)
    ax.hist(samples, bins=bins, alpha=0.7, density=True, label='Samples')
    
    # Overlay theoretical distribution
    x_theory = np.logspace(np.log10(x_min), np.log10(x_max), 100)
    # Normalized power law
    if alpha != 1:
        norm = (1-alpha)/(x_max**(1-alpha) - x_min**(1-alpha))
        p_theory = abs(norm) * x_theory**(-alpha)  # abs() for alpha < 1
    else:
        norm = 1/np.log(x_max/x_min)
        p_theory = norm / x_theory
    
    ax.plot(x_theory, p_theory, 'r-', linewidth=2, label='Theory')
    
    ax.set_xscale('log')
    ax.set_yscale('log')
    ax.set_xlabel(r'Mass $[M_☉])
    ax.set_ylabel('Probability Density')
    ax.set_title(rf'$\alpha = ${alpha}')
    ax.legend()
    ax.grid(True, alpha=0.3)
    
    # Add statistics
    mean_mass = np.mean(samples)
    median_mass = np.median(samples)
    text_str = rf'Mean: {mean_mass:.1f} $M_☉ + '\n' + rf'Median: {median_mass:.1f} $M_☉
    ax.text(0.05, 0.95, text_str, transform=ax.transAxes, 
            verticalalignment='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))

plt.suptitle('Power Law Sampling with Different Slopes', fontsize=14)
plt.tight_layout()
plt.show()

# Show how the distribution changes with alpha
print("Effect of power law slope α on stellar mass distribution:")
print(f"{'α':<5} {'Mean [M_☉]':<12} {'Median [M_☉]':<12} {'% below 1 M_☉':<15}")
print("-" * 50)

for alpha in alphas:
    samples = sample_simple_power_law(alpha, x_min, x_max, 10000)
    mean_mass = np.mean(samples)
    median_mass = np.median(samples)
    frac_below_1 = np.sum(samples < 1.0) / len(samples) * 100
    print(f"{alpha:<5.1f} {mean_mass:<12.2f} {median_mass:<12.2f} {frac_below_1:<15.1f}")

Why the Inverse Transform is Powerful¶

The inverse transform method is the foundation of Monte Carlo sampling because:

1. It’s exact: Every uniform random number produces exactly one sample from your distribution

2. It’s efficient: 100% of random numbers are used (no rejections)

3. It’s deterministic: The same uniform input always gives the same output (useful for debugging)

4. It works in any dimension: Can be extended to multivariate distributions

However, it requires that you can:

• Compute the CDF analytically (or numerically)

• Invert the CDF (analytically or numerically)

When these conditions aren’t met, you’ll need other methods like rejection sampling (Section 3.6) or Markov Chain Monte Carlo (Project 4).

Example: Exponential distribution (photon path lengths!)

• PDF: $f(x) = \lambda e^{-\lambda x}$

• CDF: $F(x) = 1 - e^{-\lambda x}$

• Inverse: $x = -\frac{1}{\lambda}\ln(1-u)$

Simple implementation:

# Sample exponential distribution (e.g., photon path lengths)
def sample_exponential(lambda_param, n_samples):
    u = np.random.uniform(0, 1, n_samples)
    # Inverse CDF transform
    return -np.log(1 - u) / lambda_param

# Example: mean free path of 1 pc
mfp = 1.0  # pc
lambda_param = 1.0 / mfp
path_lengths = sample_exponential(lambda_param, 10000)
print(f"Mean path: {np.mean(path_lengths):.2f} pc")
print(f"Theory: {mfp:.2f} pc")

4.6 Rejection Sampling¶

Priority: 🟡 Standard Path When inverse transform fails, use rejection sampling:

Example: Throw darts at rectangle, keep those under curve.

# Example: Sample from arbitrary distribution
def rejection_sample(f, x_min, x_max, f_max, n_samples):
    """Sample using rejection method"""
    samples = []
    n_tries = 0

    while len(samples) < n_samples:
        x = np.random.uniform(x_min, x_max)
        y = np.random.uniform(0, f_max)
        n_tries += 1

        if y <= f(x):
            samples.append(x)

    efficiency = n_samples / n_tries
    print(f"Efficiency: {efficiency:.1%}")
    return np.array(samples)

4.7 Spatial Distributions: The Plummer Sphere¶

Priority: 🔴 Essential for Project 2 Now that you can sample stellar masses from the Kroupa IMF, you need to place these stars in space. The Plummer sphere is a standard density profile for globular clusters and stellar systems.

The Plummer Profile¶

The Plummer sphere has density profile:

where:

• $M_{total}$ is the total cluster mass

• $a$ is the Plummer radius (scale length)

• $r$ is the distance from center

This creates a smooth, centrally concentrated distribution that avoids the infinite density at r=0 that plagues some other profiles.

The Sampling Challenge¶

You can’t sample directly from $\rho(r)$ because density isn’t probability! The key insight: the probability of finding a star between radius $r$ and $r+dr$ is proportional to the mass in that shell, not the density.

So the radial PDF is:

Approach 1: Mass Profile Method¶

The enclosed mass within radius $r$ is:

The cumulative mass fraction (which IS a CDF!) is:

To sample:

1. Generate $u \sim U(0,1)$ (this represents a mass fraction)

2. Solve $u = \frac{r^3}{(r^2 + a^2)^{3/2}}$ for $r$

3. Place the star at radius $r$ with random angular coordinates

Approach 2: Rejection Sampling in 3D¶

Alternatively, you could:

1. Generate random positions in a cube containing your cluster

2. Accept/reject based on the density at that position

3. But this becomes inefficient for large clusters!

From Radius to 3D Positions¶

Once you have radius $r$:

1. Generate random direction on sphere:

   • $\theta \sim U(0, \pi)$ (polar angle)

   • $\phi \sim U(0, 2\pi)$ (azimuthal angle)

2. Convert to Cartesian:

   • $x = r \sin\theta \cos\phi$

   • $y = r \sin\theta \sin\phi$

   • $z = r \cos\theta$

Or use the simpler method: generate a random unit vector and scale by $r$.

Verification¶

Your sampled positions should reproduce:

1. The density profile $\rho(r)$ when binned radially

2. The enclosed mass profile $M(<r)$

3. The half-mass radius $r_{half} \approx 1.3a$ for Plummer

Part 4 Synthesis: Theory to Computation¶