Probability Distribution Functions: PMF, PDF, CDF, and Common Distributions

Solution: CDF with Binary Search

class ProbabilityDistribution {
    constructor(probabilities) {
        if (!probabilities || probabilities.length === 0) {
            throw new Error('Empty probability array');
        }

        this.probabilities = probabilities;
        this.cdf = this.buildCDF(probabilities);
    }

    buildCDF(probs) {
        // Normalize if needed
        const sum = probs.reduce((a, b) => a + b, 0);
        if (Math.abs(sum - 1.0) > 1e-6) {
            probs = probs.map(p => p / sum);
        }

        // Build cumulative distribution
        const cdf = [];
        let cumulative = 0;

        for (const p of probs) {
            cumulative += p;
            cdf.push(cumulative);
        }

        // Ensure last value is exactly 1.0 (handle float errors)
        cdf[cdf.length - 1] = 1.0;

        return cdf;
    }

    sample() {
        const r = Math.random();

        // Binary search in CDF
        let left = 0;
        let right = this.cdf.length - 1;

        while (left < right) {
            const mid = Math.floor((left + right) / 2);

            if (this.cdf[mid] < r) {
                left = mid + 1;
            } else {
                right = mid;
            }
        }

        return left;
    }

    // Sample multiple times and return frequency distribution
    sampleMultiple(n) {
        const counts = Array(this.probabilities.length).fill(0);

        for (let i = 0; i < n; i++) {
            counts[this.sample()]++;
        }

        // Convert to percentages
        return counts.map(c => (c / n * 100).toFixed(2) + '%');
    }
}

// Test
const probs = [0.1, 0.3, 0.4, 0.2];
const dist = new ProbabilityDistribution(probs);

console.log('Expected probabilities:', probs.map(p => (p*100) + '%'));
console.log('Sampled 100,000 times:', dist.sampleMultiple(100000));

/* Expected output:
Expected probabilities: ['10%', '30%', '40%', '20%']
Sampled 100,000 times: ['10.01%', '29.98%', '40.02%', '19.99%']
*/

Alternative: Linear Search (Simple)

function sampleLinear(probabilities) {
    const r = Math.random();
    let cumulative = 0;

    for (let i = 0; i < probabilities.length; i++) {
        cumulative += probabilities[i];
        if (r < cumulative) {
            return i;
        }
    }

    // Should never reach here if probs sum to 1
    return probabilities.length - 1;
}

// Test
const probs = [0.1, 0.3, 0.4, 0.2];
const counts = [0, 0, 0, 0];

for (let i = 0; i < 100000; i++) {
    counts[sampleLinear(probs)]++;
}

console.log('Linear search results:');
counts.forEach((c, i) => {
    console.log(`Index ${i}: ${(c/100000*100).toFixed(2)}% (expected ${probs[i]*100}%)`);
});

Python Version with Chi-Square Test

import random
import bisect
from scipy import stats

class ProbabilityDistribution:
    def __init__(self, probabilities):
        self.probabilities = probabilities

        # Normalize if needed
        total = sum(probabilities)
        if abs(total - 1.0) > 1e-6:
            probabilities = [p / total for p in probabilities]

        # Build CDF
        self.cdf = []
        cumulative = 0
        for p in probabilities:
            cumulative += p
            self.cdf.append(cumulative)

        self.cdf[-1] = 1.0  # Fix float errors

    def sample(self):
        r = random.random()
        # bisect_left finds insertion point
        return bisect.bisect_left(self.cdf, r)

    def sample_multiple(self, n):
        counts = [0] * len(self.probabilities)
        for _ in range(n):
            counts[self.sample()] += 1
        return counts


# Test with statistical validation
probs = [0.1, 0.3, 0.4, 0.2]
dist = ProbabilityDistribution(probs)

n_samples = 100000
samples = dist.sample_multiple(n_samples)

print("Probability Distribution Sampling")
print("=" * 50)
print(f"Sample size: {n_samples}n")

print("IndextExpectedtObservedtDiff")
for i, (expected, observed) in enumerate(zip(probs, samples)):
    expected_count = expected * n_samples
    observed_pct = observed / n_samples
    diff = abs(observed - expected_count)
    print(f"{i}t{expected*100:.1f}%tt{observed_pct*100:.2f}%tt{diff:.0f}")

# Chi-square goodness of fit test
expected_counts = [p * n_samples for p in probs]
chi2, p_value = stats.chisquare(samples, expected_counts)

print(f"nChi-square statistic: {chi2:.4f}")
print(f"P-value: {p_value:.4f}")
print(f"Result: {'PASS' if p_value > 0.05 else 'FAIL'} (α=0.05)")

Complexity Analysis

When to use each:

• Linear: One-time sampling, small n
• Binary: Multiple samples, medium n (most common)
• Alias: Millions of samples, worth complex setup

Common Mistakes

1. Not normalizing: If probs sum to 0.9, last value never selected
2. Float precision: Use epsilon comparisons, not exact equality
3. Wrong binary search: Need first value >= r, not just >=
4. Not handling edge cases: Empty array, single element, etc.
5. Inefficient for repeated sampling: Rebuilding CDF each time

Follow-up Questions

1. How to sample from continuous distribution (e.g. Normal)? Use inverse CDF method or Box-Muller transform
2. Sample without replacement? Track used indices, adjust remaining probabilities
3. Weighted random with changing weights? Use heap or segment tree for O(log n) updates
4. Prove your sampling is correct? Chi-square test, KS test
5. Stream of probabilities? Use reservoir sampling

Real-World Applications

• Monte Carlo simulation: Sampling from probability models
• Game development: Loot drops, enemy spawns
• A/B testing: Assigning users to experiments
• Load balancing: Weighted server selection
• ML training: Sampling training examples

Related Problems

• Write a Function for r and 5() - Random sampling
• Coin Rolls - Expected value
• Red Marbles, Blue Marbles - Probability calculations