Verwerfungsmethode

Einfache Sprache

Durch die Verwerfungsmethode kann eine Stichprobe von einer beliebigen Verteilung gezogen werden.

Def. Verwerfungsmethode

Bei der Verwerfungsmethode geht man wie folgt vor, um eine Stichprobe der Größe $n$, aus dem Definitionsbereich $[a,b]$ mit $a,b\in\mathbb R$ und $a

Ziehe eine Zufallszahl aus der Uniformverteilung $[a,b]$.

Ziehe eine Zufallszahl aus der Uniformverteilung $0,\max f]$.

Falls $f(x) > y$ behalte die Zufallszahl, sonst verwerfe sie.

Wiederhole 1.-3. bis $n$ Zahlen behalten wurden.

Algorithmus

 1import numpy as np
 2
 3def rejection_sampling(f: callable, b: float, max_iter: int = 10000) -> float:
 4    """
 5    Performs rejection sampling to generate samples from distribution f.
 6    
 7    Args:
 8        f: Target probability density function (must be non-negative)
 9        b: Interval [-b, b] around zero to sample from
10        max_iter: Maximum number of iterations to prevent infinite loops
11        
12    Returns:
13        A sample x drawn from distribution f
14        
15    Raises:
16        RuntimeError: If fails to generate a valid sample within max_iter attempts
17        
18    The algorithm generates uniform proposals in [-b, b] and accepts them with
19    probability proportional to f(x). Requires knowing the maximum of f in [-b, b].
20    """
21    # Find maximum of f in the interval [-b, b]
22    x_vals = np.linspace(-b, b, 1000)
23    f_max = max(f(x) for x in x_vals)
24    
25    for _ in range(max_iter):
26        x = np.random.uniform(-b, b)
27        y = np.random.uniform(0, f_max)
28        
29        if y < f(x):
30            return x
31    
32    raise RuntimeError(f"Failed to generate sample after {max_iter} attempts")
33
34# Example usage:
35if __name__ == "__main__":
36    # Define a simple normal distribution-like function
37    def normal_pdf(x):
38        return np.exp(-x**2/2) / np.sqrt(2*np.pi)
39    
40    sample = rejection_sampling(normal_pdf, b=5)
41    print(f"Generated sample: {sample:.4f}")

Home: