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

4.1 Why Random Sampling Matters¶

Priority: 🔴 Essential

Sampling is how we generate distributions; moments are how we audit them. 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? Sampling feels like randomness; precisely, it is measure-preserving transforms from uniform draws. Catchphrase -> Precision: “Random draws solve deterministic problems” means unbiased estimators converge with standard error $\propto 1/\sqrt{N}$ under the right assumptions.

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

• Statistical understanding (Parts 1-2)

• Computational modeling (your projects)

{#fig-monte-carlo-pi-estimation width=“100%”}

What to notice:

• The estimate fluctuates stochastically around the true value.

• Error decreases like $1/\sqrt{N}$, not linearly with $N$.

• Accuracy gains require rapidly increasing sample counts.

Project Hook: This appears in Project 2 when you transform uniform RNG draws into physically consistent mass, radius, and velocity samples.

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} p(x')\,dx'}$

While the probability density function (pdf) $p(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 $p(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(\infty) = 1.0$ means 100% of stars have some finite mass

Key properties of CDFs:

• Always non-decreasing: $F(x_1) \le F(x_2)$ if $x_1 < x_2$

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

• Continuous from the right: Important for discrete distributions

• Derivative gives pdf: $dF/dx = p(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.

Assumption break cue: This method can fail in practice when the CDF is not reliably invertible, numerical inversion is unstable, or endpoint handling is ignored.

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: $p(x) = \lambda e^{-\lambda x}$ for x $\ge$ 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)$

Equivalent form: $x = -\frac{1}{\lambda}\ln(u)$, since $u$ and $(1-u)$ are both uniform on [0,1].

import numpy as np

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

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

Copy/paste quick-check: Imports: numpy. Seed: set rng = np.random.default_rng(536) for repeatable means. Variability: sample mean fluctuates around the theoretical mean at finite draw counts.

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: $p(x) = \frac{(\alpha-1) x^{-\alpha}}{x_{\min}^{1-\alpha} - x_{\max}^{1-\alpha}}$ (for $\alpha \neq 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 $\alpha = 1$: The integral of $1/x$ is $\ln(x)$, so: $\boxed{x = x_{\min} \left(\frac{x_{\max}}{x_{\min}}\right)^u}$

{#fig-power-law-sampling width=“100%”}

Implementation and Visualization¶

The figure above shows our power law sampling in action. Keep the production path small: sample, then run the Sampling Diagnostic Trio.

import numpy as np

def sample_simple_power_law(alpha, x_min, x_max, n_samples, rng=None):
    if rng is None:
        rng = np.random.default_rng()
    u = rng.uniform(0, 1, n_samples)
    return (x_min**(1-alpha) + u*(x_max**(1-alpha) - x_min**(1-alpha)))**(1/(1-alpha))

samples = sample_simple_power_law(2.35, 0.1, 100, 10_000, rng=np.random.default_rng(536))
print(np.mean(samples), np.var(samples), np.mean(samples > 10.0))

Use the Sampling Diagnostic Trio here: histogram-vs-theory for pdf shape, empirical CDF (or QQ) for distributional agreement, and moments plus a high-mass tail fraction for inference.

import numpy as np
import matplotlib.pyplot as plt

def sample_simple_power_law(alpha, x_min, x_max, n_samples, rng=None):
    """
    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
    """
    if rng is None:
        rng = np.random.default_rng()

    u = rng.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
rng = np.random.default_rng(536)

for ax, alpha in zip(axes, alphas):
    # Generate samples
    samples = sample_simple_power_law(alpha, x_min, x_max, 10000, rng=rng)
    
    # 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("Mass [M_sun]")
    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 = f"Mean: {mean_mass:.1f} M_sun\nMedian: {median_mass:.1f} M_sun"
    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 alpha on stellar mass distribution:")
print(f"{'alpha':<8} {'Mean [M_sun]':<14} {'Median [M_sun]':<16} {'% below 1 M_sun':<18}")
print("-" * 50)

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

Copy/paste quick-check: Imports: numpy, matplotlib. Seed: set rng = np.random.default_rng(536) for reproducible sampling tables/plots. Variability: slope-dependent trends persist, while exact summary values vary with seed and sample count.

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 4.3) or Markov Chain Monte Carlo (Project 4).

Example: Exponential distribution (photon path lengths!)

• pdf: $p(x) = \lambda e^{-\lambda x}$

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

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

• Equivalent form: $x = -\frac{1}{\lambda}\ln(u)$ because $u$ and $(1-u)$ are both Uniform$(0,1)$

Simple implementation:

import numpy as np

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

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

Copy/paste quick-check: Imports: numpy. Seed: set rng = np.random.default_rng(536) to reproduce output means. Variability: finite-sample mean fluctuates around the target mean free path.

4.3 Rejection Sampling¶

Priority: 🟡 Standard Path

When inverse transform fails, use rejection sampling:

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

import numpy as np

# Example: Sample from arbitrary distribution
def rejection_sample(p_func, x_min, x_max, f_max, n_samples, max_tries=10_000_000, rng=None):
    """Sample using rejection method"""
    if rng is None:
        rng = np.random.default_rng()
    samples = []
    n_tries = 0

    while len(samples) < n_samples:
        if n_tries >= max_tries:
            raise RuntimeError("Exceeded max_tries; check f_max envelope or target PDF.")
        x = rng.uniform(x_min, x_max)
        y = rng.uniform(0, f_max)
        n_tries += 1

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

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

Copy/paste quick-check: Imports: numpy. Seed recommendation: pass rng=np.random.default_rng(536) for repeatable envelope comparisons. Variability: efficiency varies stochastically even with fixed target PDF.

{#fig-random-sampling-methods width=“100%”}

What to notice:

• Inverse transform is exact and efficient when an invertible CDF exists.

• Rejection sampling is flexible but wastes draws depending on envelope tightness.

• The same uniform RNG can support very different target distributions.

4.4 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

Assumption break cue: Radial sampling is biased if isotropy is violated in angle draws or if finite-radius truncation is introduced without updating the CDF.

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:

   • $\mu = \cos\theta \sim U(-1,1)$

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

   • $\theta = \arccos(\mu)$

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

4. The Sampling Diagnostic Trio: histogram-vs-theory for $\rho(r)$, empirical CDF-vs-theory for $M(<r)/M_{\mathrm{total}}$, and summary moments/tail fraction for radial structure

{#fig-plummer-sphere-sampling width=“100%”}

What to notice:

• Radial sampling matches both $\rho(r)$ and $M(<r)$ diagnostics.

• The half-mass radius is recovered near the theoretical value.

• Isotropic angle sampling is required to avoid directional bias.

Part 4 Synthesis: Theory to Computation¶