def adaptive_rejection_sampling(logpdf: callable,<EOL>a: float, b: float,<EOL>domain: Tuple[float, float],<EOL>n_samples: int,<EOL>random_stream=None):
assert(hasattr(logpdf, "<STR_LIT>"))<EOL>assert(len(domain) == <NUM_LIT:2>), "<STR_LIT>"<EOL>assert(domain[<NUM_LIT:1>] >= domain[<NUM_LIT:0>]), "<STR_LIT>"<EOL>assert(n_samples >= <NUM_LIT:0>), "<STR_LIT>"<EOL>if random_stream is None:<EOL><INDENT>random_stream = RandomState()<EOL><DEDENT>if a >= b or isinf(a) or isinf(b) or a < domain[<NUM_LIT:0>] or b > domain[<NUM_LIT:1>]:<EOL><INDENT>raise ValueError("<STR_LIT>")<EOL><DEDENT>n_derivative_steps = <NUM_LIT> * (b - a)<EOL>S = (a, a + n_derivative_steps, b - n_derivative_steps, b)<EOL>if domain[<NUM_LIT:0>] == float("<STR_LIT>"):<EOL><INDENT>derivative_sign = sign(logpdf(a + n_derivative_steps) - logpdf(a))<EOL>positive_derivative = derivative_sign > <NUM_LIT:0><EOL>assert(positive_derivative), "<STR_LIT>"<EOL><DEDENT>if domain[<NUM_LIT:1>] == float("<STR_LIT>"):<EOL><INDENT>derivative_sign = sign(logpdf(b) - logpdf(b - n_derivative_steps))<EOL>negative_derivative = derivative_sign < <NUM_LIT:0><EOL>assert(negative_derivative), "<STR_LIT>"<EOL><DEDENT>n_initial_mesh_points = <NUM_LIT:3><EOL>S = unique(<EOL>(S[<NUM_LIT:0>], *(linspace(S[<NUM_LIT:1>], S[<NUM_LIT:2>], num=n_initial_mesh_points + <NUM_LIT:2>)), S[<NUM_LIT:3>])<EOL>)<EOL>fS = tuple(logpdf(s) for s in S)<EOL>lower_hull, upper_hull = compute_hulls(S=S, fS=fS, domain=domain)<EOL>samples = []<EOL>while len(samples) < n_samples:<EOL><INDENT>mesh_changed = False<EOL>x = sample_upper_hull(upper_hull, random_stream=random_stream)<EOL>lh_val, uh_val = evaluate_hulls(x, lower_hull, upper_hull)<EOL>U = random_stream.rand()<EOL>if log(U) <= lh_val - uh_val:<EOL><INDENT>samples.append(x)<EOL><DEDENT>elif log(U) <= logpdf(x) - uh_val:<EOL><INDENT>samples.append(x)<EOL>mesh_changed = True<EOL><DEDENT>else:<EOL><INDENT>mesh_changed = True<EOL><DEDENT>if mesh_changed:<EOL><INDENT>S = sorted([*S, x])<EOL>fS = tuple(logpdf(s) for s in S)<EOL>lower_hull, upper_hull = compute_hulls(S=S, fS=fS, domain=domain)<EOL><DEDENT><DEDENT>return samples<EOL>
Adaptive rejection sampling samples exactly (all samples are i.i.d) and efficiently from any univariate log-concave distribution. The basic idea is to successively determine an envelope of straight-line segments to construct an increasingly accurate approximation of the logarithm.
It does not require any normalization of the target distribution.

Parameters
----------
logpdf: callable
    Univariate function that computes :math:`log(f(u))`
    for a given :math:`u`, where :math:`f(u)` is proportional
    to the target density to sample from.

a: float
    Lower starting point used to initialize the hulls.
    Must lie in the domain of the logpdf and it
    must hold: :math:`a < b`.

b: float
    Upper starting point used to initialize the hulls.
    Must lie in the domain of the logpdf and it
    must hold: :math:`a < b`.

domain : Tuple[float, float]
    Domain of `logpdf`.
    May be unbounded on either or both sides,
    in which case `(float("-inf"), float("inf"))`
    would be passed.
    If this domain is unbounded to the left,
    the derivative of the logpdf
    for x<= a must be positive.
    If this domain is unbounded to the right the derivative of the logpdf for x>=b
    must be negative.

n_samples: int
    Number of samples to draw.

random_stream : RandomState, optional
    Seeded random number generator object with same interface as a NumPy
    RandomState object. Defaults to `None` in which case a NumPy
    RandomState seeded from `/dev/urandom` if available or the clock if not
    will be used.

Returns
----------
samples : list
    A list of samples drawn from the
    target distribution :math:`f`
    with the given `logpdf`.

Examples
----------
Sampling from a simple gaussian, adaptive rejection sampling style.
We use the logpdf of a standard gaussian and this small code snippet
demonstrates that our sample approximation accurately approximates the mean:

>>> from math import isclose
>>> from numpy import log, exp, mean
>>> gaussian_logpdf = lambda x, sigma=1: log(exp(-x ** 2 / sigma))
>>> a, b = -2, 2  # a < b must hold
>>> domain = (float("-inf"), float("inf"))
>>> n_samples = 10000
>>> samples = adaptive_rejection_sampling(logpdf=gaussian_logpdf, a=a, b=b, domain=domain, n_samples=n_samples)
>>> isclose(mean(samples), 0.0, abs_tol=1e-02)
True