How Monte Carlo Method Works: Implementation Guide

Introduction and Core Concepts

The Monte Carlo method is a broad class of computational algorithms that rely on repeated random sampling to obtain numerical results. Named after the famous Monte Carlo Casino in Monaco due to its inherent reliance on chance, the technique was formally developed in the 1940s by scientists working on the Manhattan Project, including Stanislaw Ulam and John von Neumann. Faced with intractable problems in neutron diffusion that could not be solved analytically, they realized that random simulations could provide approximate yet remarkably accurate solutions.

At its core, the Monte Carlo method transforms a deterministic problem into a probabilistic one. Instead of trying to solve complex integrals, differential equations, or optimization problems directly, we simulate a large number of random scenarios and aggregate the results. The fundamental principle underpinning this approach is the Law of Large Numbers from probability theory, which states that as the number of independent trials increases, the average of the results obtained will converge to the expected value. Similarly, the Central Limit Theorem allows us to quantify the error and confidence intervals of our estimates.

Why does the Monte Carlo method matter today? In an era of massive computational power, many real-world problems remain analytically unsolvable or computationally prohibitive for deterministic methods. High-dimensional integration, complex risk modeling, and systems with inherent uncertainty are prime examples. Monte Carlo methods provide a flexible, scalable, and often surprisingly simple alternative. They are embarrassingly parallel, meaning they can be easily distributed across thousands of cores or machines, making them ideal for modern cloud and high-performance computing environments. From estimating the value of pi to pricing exotic financial derivatives and training artificial intelligence, the method's universality makes it an essential tool for any engineer, data scientist, or quantitative analyst.

The Core Philosophy: Randomness as a Computational Tool

Traditional algorithms follow a precise, deterministic path from input to output. Monte Carlo algorithms embrace randomness as a feature, not a bug. By sampling from a probability distribution that describes the problem domain, we can explore the solution space efficiently without having to enumerate every possibility. This is particularly powerful when the solution space is vast or continuous.

A key concept is the distinction between Monte Carlo and Las Vegas algorithms. A Monte Carlo algorithm may produce an incorrect result with a small probability, but its runtime is typically fixed or predictable. In contrast, a Las Vegas algorithm always produces a correct result, but its runtime is probabilistic. Most practical implementations of the Monte Carlo method fall into the first category, where we trade a tiny, quantifiable error probability for significant gains in simplicity and speed.

How It Works

The mechanism of the Monte Carlo method can be broken down into a general four-step workflow that applies to nearly every application. First, we define a domain of possible inputs. Second, we generate random inputs from a probability distribution over that domain. Third, we perform a deterministic computation on each input. Finally, we aggregate the results to approximate the desired quantity. This simple pattern is deceptively powerful.

Consider the classic example of estimating the value of pi. We know that a circle of radius r inscribed in a square with side length 2r has an area ratio of pi/4. If we randomly throw darts at the square, the proportion of darts that land inside the circle will converge to pi/4. By generating thousands or millions of random (x, y) points uniformly distributed in the square and counting how many satisfy x^2 + y^2 <= r^2, we can estimate pi as 4 * (points_inside / total_points). No geometry beyond the distance formula is required.

The accuracy of a Monte Carlo estimate is governed by the number of samples, N. The standard error of the mean decreases at a rate of 1 / sqrt(N), regardless of the dimensionality of the problem. This is a critical advantage. While deterministic numerical integration methods like Simpson's rule suffer from the curse of dimensionality — their error grows exponentially with dimensions — Monte Carlo's convergence rate is dimension-independent. To halve the error, you need to quadruple the number of samples. This O(1/sqrt(N)) convergence is slow but predictable and robust.

Probability Distributions and Sampling

The choice of sampling distribution is crucial. The simplest case is uniform sampling, where every point in the domain is equally likely. However, many problems benefit from importance sampling, where we sample more frequently from regions that contribute more to the final result. For example, when integrating a function that is sharply peaked in one area, uniform sampling wastes many samples in low-value regions. By sampling from a distribution that mimics the shape of the function, we can dramatically reduce variance.

Generating high-quality pseudo-random numbers is the engine of any Monte Carlo simulation. Modern programming languages use algorithms like Mersenne Twister or PCG to produce sequences that are statistically indistinguishable from true randomness for most purposes. The quality of the random number generator directly impacts the quality of the result, especially for simulations requiring billions of samples. For cryptographic or ultra-high-precision scientific work, hardware-based true random number generators may be used.

Convergence and Error Estimation

Unlike deterministic algorithms that produce an exact answer, Monte Carlo methods produce an estimate with an associated confidence interval. Using the sample variance, we can compute the standard error and construct intervals such as estimate ± 1.96 * standard_error for a 95% confidence level. This built-in error quantification is a major strength, as it tells us how much we can trust our result and how many more samples are needed to reach a desired precision.

Implementation

Implementing a Monte Carlo method in Python is straightforward thanks to the powerful random and numpy libraries. NumPy is strongly preferred for performance because it allows vectorized operations on large arrays, avoiding slow Python loops. The general structure involves setting a random seed for reproducibility, generating a large array of random samples, applying a vectorized function, and averaging the results.

Below is a foundational implementation for estimating pi. Note how NumPy generates all random points at once and uses boolean masking to count points inside the circle efficiently. This vectorized approach can be millions of times faster than a naive Python loop for large N.

import numpy as np
import math

def estimate_pi_vectorized(num_samples: int, seed: int = 42) -> float:
    """
    Estimate pi using Monte Carlo method with vectorized NumPy operations.
    
    Args:
        num_samples: Number of random points to generate
        seed: Random seed for reproducibility
    
    Returns:
        Estimated value of pi
    """
    # Set seed for reproducible results
    np.random.seed(seed)
    
    # Generate random points uniformly in [-1, 1] x [-1, 1]
    # Shape: (num_samples, 2)
    points = np.random.uniform(-1, 1, size=(num_samples, 2))
    
    # Calculate distance from origin for each point
    # Using vectorized operation: sqrt(x^2 + y^2)
    distances_squared = np.sum(points**2, axis=1)
    
    # Count points inside unit circle (distance < 1)
    points_inside = np.sum(distances_squared < 1.0)
    
    # Area ratio: circle_area / square_area = pi / 4
    # So pi ~= 4 * (points_inside / num_samples)
    pi_estimate = 4 * points_inside / num_samples
    
    return pi_estimate

# Example usage
for n in [1000, 10000, 100000, 1000000]:
    estimate = estimate_pi_vectorized(n)
    error = abs(estimate - math.pi)
    print(f"N={n:7d} | Estimate={estimate:.6f} | Error={error:.6f}")

A more general application is Monte Carlo integration. Any definite integral can be interpreted as an expected value. To integrate f(x) over [a, b], we can sample x uniformly from [a, b] and compute the average value of f(x) multiplied by the interval width (b - a). This works for functions that have no closed-form antiderivative.

import numpy as np
from typing import Callable

def monte_carlo_integration(
    func: Callable[[np.ndarray], np.ndarray],
    a: float,
    b: float,
    num_samples: int = 100000,
    seed: int = 0
) -> tuple[float, float]:
    """
    Generic Monte Carlo integrator for 1D functions.
    
    Estimates integral of func from a to b.
    Returns (estimate, standard_error).
    """
    np.random.seed(seed)
    
    # Generate uniform samples in [a, b]
    samples = np.random.uniform(a, b, size=num_samples)
    
    # Evaluate function at all sample points (vectorized)
    func_values = func(samples)
    
    # Monte Carlo estimate: (b - a) * mean(f(x))
    interval_width = b - a
    estimate = interval_width * np.mean(func_values)
    
    # Estimate standard error for confidence interval
    # std_error = (b - a) * std(f(x)) / sqrt(N)
    std_dev = np.std(func_values, ddof=1)
    std_error = interval_width * std_dev / np.sqrt(num_samples)
    
    return estimate, std_error

# Example: Integrate f(x) = x^2 from 0 to 2
# Analytical result is 8/3 ~= 2.6667
def f(x: np.ndarray) -> np.ndarray:
    return x**2

estimate, error = monte_carlo_integration(f, 0, 2, num_samples=500000)
print(f"Estimate: {estimate:.5f} +/- {1.96*error:.5f} (95% CI)")
print(f"True value: {8/3:.5f}")

# Example: Integrate sin(x) from 0 to pi (true value = 2)
def g(x: np.ndarray) -> np.ndarray:
    return np.sin(x)

est2, err2 = monte_carlo_integration(g, 0, np.pi, num_samples=500000)
print(f"Integral of sin(x) from 0 to pi: {est2:.5f} +/- {1.96*err2:.5f}")

For robust production code, always consider reproducibility, error reporting, and performance. Setting a seed allows others to reproduce your exact results. Returning both the estimate and its standard error provides crucial context about reliability. Using NumPy's vectorization and avoiding explicit Python loops with for i in range(N): where i < N is checked each iteration will ensure your simulation scales to millions or billions of samples.

Step-by-Step Example

Let us walk through a complete, concrete example: estimating the integral of f(x) = exp(-x^2) from 0 to 1. This is a part of the Gaussian integral and has no elementary antiderivative, making it a perfect candidate for Monte Carlo. The true value is approximately 0.746824. We will use N = 100,000 samples and track the convergence.

Step 1: Define the problem domain. Our interval is [0, 1], so the width is 1.0. This simplifies our estimate to just the mean of f(x) over the samples. Step 2: Generate random samples. We draw 100,000 values uniformly from [0, 1]. Each sample represents a random x-coordinate. Step 3: Evaluate the function. For each x, we compute exp(-x^2). This gives us the height of the function at that random point. Step 4: Aggregate. We average all 100,000 heights. Because the interval width is 1, this average is our integral estimate.

To understand convergence, we can compute a running estimate. This shows how the estimate stabilizes as N grows. Early on with few samples, the estimate fluctuates wildly. As N increases into the tens of thousands, the fluctuations dampen and the estimate settles near the true value. This behavior directly illustrates the 1/sqrt(N) error reduction.

import numpy as np
import math

def step_by_step_gaussian_integral(num_samples: int = 100000):
    """
    Walk through Monte Carlo integration of exp(-x^2) from 0 to 1
    with convergence tracking.
    """
    np.random.seed(123)
    
    # True value for comparison (using math.erf)
    # integral exp(-x^2) dx from 0 to 1 = sqrt(pi)/2 * erf(1)
    true_value = math.sqrt(math.pi) / 2 * math.erf(1)
    print(f"True value: {true_value:.6f}\n")
    
    # Step 2: Generate all samples at once
    samples = np.random.uniform(0, 1, size=num_samples)
    
    # Step 3: Evaluate function
    values = np.exp(-samples**2)
    
    # Step 4: Track convergence at checkpoints
    checkpoints = [100, 1000, 10000, 50000, 100000]
    print(f"{'N':>8} | {'Estimate':>10} | {'Error':>10} | {'Std Error':>10}")
    print("-" * 50)
    
    for n in checkpoints:
        if n > num_samples:
            continue
        # Running mean of first n values
        current_estimate = np.mean(values[:n])
        current_std = np.std(values[:n], ddof=1)
        std_error = current_std / math.sqrt(n)
        error = abs(current_estimate - true_value)
        print(f"{n:8d} | {current_estimate:10.6f} | {error:10.6f} | {std_error:10.6f}")

    # Final 95% confidence interval
    final_estimate = np.mean(values)
    final_std_error = np.std(values, ddof=1) / math.sqrt(num_samples)
    print(f"\nFinal 95% CI: [{final_estimate - 1.96*final_std_error:.6f}, "
          f"{final_estimate + 1.96*final_std_error:.6f}]")
    print(f"Interval contains true value? "
          f"{(final_estimate - 1.96*final_std_error < true_value < final_estimate + 1.96*final_std_error)}")

step_by_step_gaussian_integral(100000)

Running this code typically produces output where the error drops from ~0.02 at N=100 to ~0.001 at N=100,000. The standard error column shows our internal estimate of uncertainty shrinking predictably. Importantly, the final confidence interval almost always contains the true value, validating our error quantification. If you need higher precision, you know exactly how many more samples are required: to reduce error by a factor of 10, you need 100 times more samples.

This example also highlights a key debugging practice: always test your Monte Carlo implementation on a problem with a known analytical solution before applying it to an unknown problem. By comparing against the known value of erf(1), we verify our sampler, function evaluation, and aggregation logic are all correct.

Complexity Analysis

Analyzing Monte Carlo methods requires a shift in perspective from traditional deterministic complexity. The runtime is almost always linear in the number of samples, but the more interesting question is how many samples are needed to achieve a target accuracy, epsilon. This sample complexity is what determines the practical cost.

Time complexity is straightforward: generating N samples and evaluating the function N times takes O(N) time. Each sample is independent, so the work scales linearly. Space complexity is O(1) if done iteratively or O(N) if all samples are stored at once for vectorization. In practice, we often use vectorized batches to balance memory and speed, processing say 1 million samples at a time in a loop.

The error complexity is O(1 / sqrt(N)) for standard Monte Carlo. To achieve an absolute error of epsilon, we need N = O(1 / epsilon^2) samples. This is independent of dimension d, which is the method's superpower. Deterministic grid-based integration with error epsilon requires N = O(1 / epsilon^d) points, which becomes astronomically large for d > 4. This is why Monte Carlo is the only feasible approach for high-dimensional integrals in physics and finance.

Here, C is the cost of one function evaluation. If C is expensive, variance reduction techniques become critical to reduce N. The best, average, and worst cases are identical in terms of Big O because the algorithm's structure does not depend on input values — it always performs N iterations. The only variation is in the constant factor due to random fluctuations in convergence, but asymptotically it remains the same.

For parallel execution, the time complexity becomes O(N / P) with P processors, plus a small O(log P) aggregation overhead. This near-perfect scalability is why Monte Carlo simulations are ideal for MapReduce, Spark, and GPU computing. You can achieve linear speedup simply by distributing samples across cores.

Variations and Optimizations

Plain Monte Carlo is powerful but can be slow to converge. A suite of variance reduction techniques can dramatically improve efficiency, often reducing the number of samples needed by orders of magnitude without increasing N. These optimizations are essential for production systems.

Importance Sampling is the most widely used optimization. Instead of sampling uniformly, we sample from a proposal distribution q(x) that is large where |f(x)| is large. We then weight each sample by p(x)/q(x) to keep the estimate unbiased. If q(x) closely mimics f(x), variance can approach zero. Antithetic Variates is another simple yet effective technique: for every random sample U, we also use 1 - U. This induces negative correlation between pairs of samples, reducing overall variance, especially for monotonic functions.

Quasi-Monte Carlo (QMC) replaces pseudo-random numbers with low-discrepancy sequences like Sobol or Halton sequences. These sequences are designed to fill space more uniformly than random points, avoiding clumps and gaps. For integration, QMC can achieve convergence close to O(1/N) instead of O(1/sqrt(N)), a massive improvement. However, error estimation is more difficult because samples are no longer independent.

import numpy as np
import math

def importance_sampling_example(num_samples: int = 100000):
    """
    Compare plain Monte Carlo vs Importance Sampling
    for integrating exp(-x^2) on [0, 1] is not ideal,
    so we use a peaked function: f(x) = exp(-50*(x-0.3)^2) on [0,1]
    Plain MC wastes samples; Importance Sampling focuses on peak.
    """
    np.random.seed(0)
    
    def f(x):
        return np.exp(-50 * (x - 0.3)**2)
    
    # Plain Monte Carlo (uniform)
    uniform_samples = np.random.uniform(0, 1, size=num_samples)
    plain_values = f(uniform_samples)
    plain_est = np.mean(plain_values) * 1.0
    plain_var = np.var(plain_values) / num_samples
    
    # Importance Sampling: sample from Normal(0.3, 0.1) truncated to [0,1]
    # For simplicity, sample from normal and reject outside [0,1]
    # Proposal q(x) ~ Normal(0.3, 0.15)
    proposal_samples = np.random.normal(0.3, 0.15, size=num_samples*2)
    proposal_samples = proposal_samples[(proposal_samples >= 0) & (proposal_samples <= 1)][:num_samples]
    
    # Need to handle if we filtered too many; refill if needed
    while len(proposal_samples) < num_samples:
        extra = np.random.normal(0.3, 0.15, size=num_samples)
        extra = extra[(extra >= 0) & (extra <= 1)]
        proposal_samples = np.concatenate([proposal_samples, extra])[:num_samples]
    
    # Compute importance weights: p(x)/q(x)
    # p(x) = 1 on [0,1] (uniform target), q(x) = normal pdf
    def normal_pdf(x, mean=0.3, sigma=0.15):
        return (1 / (sigma * math.sqrt(2*math.pi))) * np.exp(-0.5 * ((x-mean)/sigma)**2)
    
    q_vals = normal_pdf(proposal_samples)
    p_vals = np.ones_like(proposal_samples)  # uniform density = 1
    weights = p_vals / q_vals
    # Correct for truncation: approximate normalization
    # For demo, simple weight adjustment
    weighted_values = f(proposal_samples) * weights
    # Normalize weights to account for truncated proposal
    is_est = np.mean(weighted_values) * 1.0  # interval width
    # Rough variance estimate
    is_var = np.var(weighted_values) / num_samples
    
    print(f"Plain MC Estimate: {plain_est:.6f} Var: {plain_var:.8f}")
    print(f"Importance Sampling Estimate: {is_est:.6f} Var: {is_var:.8f}")
    print(f"Variance Reduction: {plain_var / is_var:.2f}x")

importance_sampling_example(200000)

Other important variants include Markov Chain Monte Carlo (MCMC) for sampling from complex, high-dimensional probability distributions where direct sampling is impossible. MCMC constructs a Markov chain whose stationary distribution is the target distribution, with Metropolis-Hastings and Gibbs sampling as the most famous algorithms. Stratified Sampling divides the domain into strata and samples within each stratum, ensuring coverage. Control Variates uses a known, similar function with a known integral to reduce variance of the unknown function.

Choosing the right variant depends on problem structure. If you know where the function is large, use Importance Sampling. If you need high accuracy in low to moderate dimensions (d < 10), try Quasi-Monte Carlo. If you need to sample from a Bayesian posterior, MCMC is the standard. Often, combining techniques — such as importance sampling with antithetic variates — yields the best results.

Real-World Applications

Monte Carlo methods are ubiquitous in industry and science precisely because they handle complexity and uncertainty so naturally. Their ability to model stochastic systems makes them indispensable where deterministic models fail.

In finance, Monte Carlo simulation is the workhorse for derivatives pricing and risk management. The price of an option depends on the future distribution of the underlying asset price, often modeled as geometric Brownian motion. By simulating thousands of price paths and averaging the discounted payoff, banks price exotic options that have no closed-form formula. Value at Risk (VaR) and Expected Shortfall are also computed via Monte Carlo by simulating portfolio losses under random market scenarios. Web-scale companies use similar simulations for A/B test power analysis and revenue forecasting under uncertainty.

In computer systems, Monte Carlo methods help with performance and reliability analysis. Database query optimizers use sampling to estimate selectivity and cardinality without scanning entire tables. Content Delivery Networks and load balancers simulate request traffic to test capacity planning and failover strategies. Network reliability analysis — estimating the probability that a network remains connected when links fail randomly — is a classic Monte Carlo application, as exact computation is #P-complete.

In artificial intelligence, Monte Carlo Tree Search (MCTS) revolutionized game playing, powering AlphaGo. MCTS randomly simulates many possible future game playouts to evaluate the promise of a move without needing a heuristic evaluation function. In Bayesian machine learning, MCMC methods like Hamiltonian Monte Carlo are used to sample from posterior distributions over model parameters, enabling uncertainty quantification.

Perhaps the most visually striking application is in computer graphics. The rendering equation that describes how light bounces in a scene is an intractable high-dimensional integral. Path tracing, the algorithm behind photorealistic images in films, uses Monte Carlo to randomly sample light paths. More samples per pixel reduce noise, which is why rendered images converge from grainy to clean as sampling progresses.

Comparison with Alternatives

Monte Carlo is not always the best tool. Understanding its trade-offs against deterministic alternatives is key to choosing the right algorithm for the job. The decision hinges on dimensionality, required accuracy, smoothness of the function, and computational budget.

Deterministic numerical integration methods like the Trapezoidal Rule, Simpson's Rule, and Gaussian Quadrature are far more efficient in one dimension. For a smooth 1D function, Simpson's rule achieves O(1/N^4) convergence, vastly superior to Monte Carlo's O(1/sqrt(N)). However, these methods require a grid. In d dimensions, they need N^d points, making them infeasible for d > 4. Monte Carlo's dimension-independent convergence makes it the winner in high dimensions.

Analytical or closed-form solutions, when available, are always preferable. They provide exact answers in O(1) time. The problem is that most real-world integrals and stochastic systems have no closed form. Similarly, dynamic programming and exact enumeration work for small discrete problems but explode combinatorially. Monte Carlo provides an approximate answer where exact methods are impossible.

Another alternative is Quasi-Monte Carlo, which as discussed improves convergence to near O(1/N) by using deterministic low-discrepancy sequences. It outperforms plain Monte Carlo for moderate dimensions and smooth functions but loses its advantage in very high dimensions or for non-smooth functions, and it complicates error estimation. For many production systems, a hybrid approach is used: start with Quasi-Monte Carlo for speed, then switch to plain Monte Carlo to get reliable confidence intervals.

Common Mistakes and Best Practices

Even though the Monte Carlo method is conceptually simple, subtle implementation errors can lead to biased results, underestimated errors, or poor performance. Awareness of these pitfalls is what separates a naive simulation from a reliable one.

A frequent mistake is using a poor random number generator or reusing the same seed unintentionally. Using Python's built-in random module in a tight loop is slow and has a shorter period than numpy.random. Failing to set a seed makes results irreproducible, while setting the same seed inside a loop (instead of once at the start) produces identical samples every iteration, completely invalidating the simulation. Another common error is integer division or incorrect scaling — for example, forgetting to multiply by (b - a) in integration or using points_inside / num_samples without the factor of 4 when estimating pi.

Underestimating variance and misreporting confidence is a serious statistical pitfall. Beginners often run a simulation once with N=10,000, report the point estimate, and ignore the standard error. Always report the confidence interval. Also, beware of pseudo-convergence: the running average may appear to stabilize but still be far from the true value if variance is high. Run multiple independent replications with different seeds to verify stability. Do not confuse the standard deviation of the samples with the standard error of the mean; the latter is what determines the accuracy of your estimate and shrinks with sqrt(N).

Best practices include: vectorize with NumPy to achieve speed; use batch processing for huge N to control memory; apply variance reduction when function evaluations are expensive; validate against known cases; and parallelize correctly by giving each worker a different, non-overlapping random seed stream (using SeedSequence in NumPy). For critical simulations, use common random numbers when comparing two system configurations — using the same random samples for both reduces variance in the difference. Finally, always visualize convergence with a running mean plot to diagnose problems early.

Checklist for Production-Ready Code

• Set seed once at the top for reproducibility, allow seed as parameter.
• Use np.random.default_rng(seed) (new Generator API) instead of legacy np.random.seed for better statistical quality.
• Return both estimate and standard error / confidence interval.
• Handle edge cases: N <= 0, invalid interval where a >= b, NaN in function values.
• Log N, estimate, error, and runtime for monitoring.
• Test with analytical cases and test that error ~ 1/sqrt(N) empirically.

Interview Questions

Monte Carlo methods are a favorite topic in interviews for quantitative, data science, and software engineering roles because they test probability, coding, and systems thinking simultaneously. Here are common questions and how to approach them.

1. How would you estimate pi using Monte Carlo? What is the error as a function of N?

Answer: Describe the dartboard method — sample (x,y) uniformly in [0,1]x[0,1] or [-1,1]x[-1,1], count fraction inside circle, multiply by 4. The estimator is p_hat = 4 * k/N where k ~ Binomial(N, pi/4). Variance is 16 * p(1-p)/N, so standard error is O(1/sqrt(N)). To get d decimal digits (~10^-d error), need N ~ 10^(2d) samples. Mention vectorized NumPy implementation and that you would report a confidence interval, not just a point estimate.

2. When would you prefer Monte Carlo integration over Simpson's rule?

Answer: In high dimensions. Simpson's rule error is O(1/N^4) in 1D but requires a grid, so in d dimensions it needs N^d points and error degrades to O(1/N^(4/d)). Monte Carlo error is O(1/sqrt(N)) independent of d. So for d >= 4-5, Monte Carlo wins. Also mention Monte Carlo needs no smoothness assumptions, handles complex domains and discontinuities better, and is trivially parallelizable, while Simpson's needs smooth derivatives and struggles with irregular boundaries.

3. Explain variance reduction. Give an example of importance sampling.

Answer: Variance reduction techniques reduce the variance of the estimator without increasing N, thus needing fewer samples for same accuracy. Importance sampling samples from a proposal distribution q(x) that concentrates samples where f(x) is large, weighting by p(x)/q(x). Example: integrating a function peaked at 0.3, uniform sampling wastes samples; sampling from Normal(0.3, sigma) focuses effort. Other techniques: antithetic variates (use U and 1-U), control variates (use known integral of similar function), stratified sampling. Emphasize that a good q(x) mimics |f(x)|.

4. How do you generate a random sample from a custom distribution, e.g., exponential, using only uniform samples?

Answer: Use inverse transform sampling. If U ~ Uniform(0,1) and F is the CDF of target distribution, then X = F^{-1}(U) has distribution F. For exponential with rate lambda, CDF is 1 - exp(-lambda x), so inverse is X = -ln(1-U)/lambda, and since 1-U is also uniform, X = -ln(U)/lambda. Mention that this requires invertible CDF; otherwise use rejection sampling or MCMC. Provide Python code: np.random.exponential internally does this, but you can implement as -np.log(np.random.uniform(0,1, size=N))/lam with check that U > 0.

5. You run a Monte Carlo simulation twice with same N and get different answers. Is your code buggy?

Answer: No, this is expected due to randomness. The estimates are random variables with variance O(1/N). They should be within a few standard errors of each other. If difference is much larger than ~3*std_error, then suspect bug, bad RNG, or insufficient N. Best practice is to report mean and confidence interval, run multiple independent replications, and check that true value (if known) lies within interval ~95% of the time. To make results reproducible, set a random seed.