diff --git a/HISTORY.rst b/HISTORY.rst
index cf47536..c292723 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -8,6 +8,8 @@ unreleased
 Added:
 ^^^^^^
 - Added functions to calculat the average total number of reflections in a room and their density. (PR #174)
+- Added a function to sample times of arrivals of reflections based on a Poisson process (PR #176)
+- Added a function to sample the amplitude or sign of reflections (PR #176)
 
 1.0.0 (2026-03-19)
 ------------------
diff --git a/pyproject.toml b/pyproject.toml
index 7aaa8b1..faffe3e 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -26,7 +26,7 @@ classifiers = [
 dependencies = [
     'pyfar>=0.8.0',
     'numpy>=1.14.0',
-    'scipy>=1.5.0',
+    'scipy>=1.6.0',
     'matplotlib',
 ]
 
diff --git a/pyrato/parametric.py b/pyrato/parametric.py
index 867ed03..dd4e485 100644
--- a/pyrato/parametric.py
+++ b/pyrato/parametric.py
@@ -4,8 +4,9 @@
 such as Sabine's theory of sound in rooms.
 """
 import numpy as np
-from typing import Union, List
+from typing import Literal, Union, List
 import pyfar as pf
+from scipy.integrate import cumulative_trapezoid
 
 
 def energy_decay_curve(
@@ -441,3 +442,274 @@ def average_number_of_reflections(
     number_of_reflections = density.time * times / 3
     return pf.TimeData(number_of_reflections, times)
 
+
+def _start_time_of_arrival_poisson_process(
+        volume : float,
+        speed_of_sound: float | None = None,
+    ) -> float:
+    """
+    The earliest time of arrival approximated as a Poisson process.
+
+    Calculated according to [#]_.
+
+    Parameters
+    ----------
+    volume : float
+        Volume of the room in m³.
+    speed_of_sound : float, None, optional
+        Speed of sound in the room. By default,
+        the :py:attr:`~pyfar.constants.reference_speed_of_sound` is used
+        which corresponds to the speed of sound in air at 20 °C.
+
+    Returns
+    -------
+    float
+        Earliest expected time of arrival in seconds.
+
+    References
+    ----------
+    .. [#] D. Schröder, “Physically based real-time auralization of
+           interactive virtual environments,” PhD Thesis, Logos-Verlag,
+           Berlin, 2011. [Online].
+           Available: https://publications.rwth-aachen.de/record/50580
+
+    """
+
+    if speed_of_sound is None:
+        speed_of_sound = pf.constants.reference_speed_of_sound
+    if speed_of_sound <= 0:
+        raise ValueError("speed_of_sound must be positive.")
+
+    if volume <= 0:
+        raise ValueError("'volume' must be positive.")
+
+    return (2*volume*np.log(2)/ (4*np.pi*speed_of_sound**3))**(1/3)
+
+
+def time_of_arrival_poisson_process(
+        volume: float,
+        times: np.ndarray,
+        speed_of_sound: float | None = None,
+        reflection_rate_limit: float = np.inf,
+        seed: int | None = None,
+    ) -> np.ndarray:
+    """Generate a time of arrival sequence based on a Poisson process.
+
+    The reflection rate is calculated using the average reflection density
+    in a diffuse sound field.
+    Note that the reflection rate increases with time, yielding a
+    non-homogeneous Poisson process. Optionally, the reflection rate can
+    be limited using the ``reflection_rate_limit`` parameter.
+    In [#]_, a maximum of 10000 reflections per second is suggested.
+
+    The implementation of the non-homogeneous Poisson process is based on
+    the transform method described in chap 5 of [#]_.
+
+    Parameters
+    ----------
+    volume : float
+        Volume of the room in m³.
+    times : numpy.ndarray
+        Time vector in seconds.
+    speed_of_sound : float, None, optional
+        Speed of sound in the room. By default, the
+        :py:attr:`~pyfar.constants.reference_speed_of_sound` is used.
+    reflection_rate_limit : float, optional
+        Maximum reflection rate in 1/s. If ``np.inf``, no limit is applied.
+        Default is ``np.inf``.
+    seed : int, None, optional
+        Seed for the random number generator. If None, a random seed is used.
+        Default is None.
+
+    Returns
+    -------
+    numpy.ndarray
+        Array of arrival times in seconds.
+
+    Examples
+    --------
+    Simulate the time of arrival of reflections in a room with a volume
+    of 100 m³ and compare the cumulative histogram to the model prediction.
+
+    .. plot::
+
+        >>> import pyrato
+        >>> import numpy as np
+        >>> import pyfar as pf
+        >>> import matplotlib.pyplot as plt
+        ...
+        >>> volume = 100
+        >>> times = np.linspace(0, .5, 200)
+        >>> toa = pyrato.parametric.time_of_arrival_poisson_process(
+        ...     volume, times)
+        ...
+        >>> plt.figure(figsize=(8, 4))
+        >>> plt.hist(
+        ...     toa, density=False, bins=50, cumulative=True, histtype='step',
+        ...     linewidth=1.5, color='C0', label='Simulation')
+        >>> ax = pf.plot.time(
+        ...     pyrato.parametric.average_number_of_reflections(
+        ...         volume, times),
+        ...     label='Model', linestyle='--',
+        ...     color='grey', linewidth=1.5)
+        >>> ax.set_ylabel('Number of reflections')
+        >>> ax.set_yscale('log')
+        >>> ax.legend(loc='lower right')
+        >>> ax.grid()
+
+    References
+    ----------
+    .. [#] D. Schröder, “Physically based real-time auralization of
+           interactive virtual environments,” PhD Thesis, Logos-Verlag,
+           Berlin, 2011. [Online].
+           Available: https://publications.rwth-aachen.de/record/50580
+    .. [#] S. M. Ross, Simulation, Sixth edition. London, United Kingdom:
+           Academic Press, 2023.
+
+    """
+
+    if speed_of_sound is None:
+        speed_of_sound = pf.constants.reference_speed_of_sound
+
+    if speed_of_sound <= 0:
+        raise ValueError("speed_of_sound must be positive.")
+
+    if volume <= 0:
+        raise ValueError("'volume' must be positive.")
+
+    if (
+        not np.isinf(reflection_rate_limit) and np.isnan(reflection_rate_limit)
+    ) or reflection_rate_limit < 0:
+        raise ValueError(
+            "'reflection_rate_limit' must be non-negative and not NaN.",
+        )
+    rng = np.random.default_rng(seed=seed)
+
+    reflection_density = average_reflection_density(
+        volume, times, speed_of_sound,
+    )
+
+    mu_values = np.minimum(
+        np.squeeze(reflection_density.time),
+        reflection_rate_limit,
+    )
+
+    mu_times = reflection_density.times
+    if np.any(np.diff(mu_times) <= 0):
+        raise ValueError("'times' must be strictly increasing.")
+    t_start = _start_time_of_arrival_poisson_process(volume, speed_of_sound)
+
+    # Cumulative intensity F(t) via numerical integration
+    cumulative_intensity = cumulative_trapezoid(mu_values, mu_times, initial=0)
+
+    # Interpolate the cumulative intensity to find the warped time values
+    F_start = np.interp(t_start, mu_times, cumulative_intensity)
+    F_end = cumulative_intensity[-1]
+
+    # expected number of arrivals
+    total_events = F_end - F_start
+
+    # Draw the total count, then place events uniformly in warped time
+    n_events = rng.poisson(total_events)
+    warped = rng.uniform(F_start, F_end, size=n_events)
+    warped.sort()
+
+    # Invert the warped time to get the arrival times by interpolation
+    arrivals = np.interp(warped, cumulative_intensity, mu_times)
+
+    return arrivals[arrivals >= t_start]
+
+
+def random_reflection_sequence(
+        arrivals : np.ndarray,
+        n_samples : int,
+        sampling_rate : float,
+        distribution : Literal['uniform', 'binary', 'normal'] = "normal",
+        seed : int | None = None,
+    ) -> pf.Signal:
+    r"""Generate reflection sequence from arrival times with random amplitudes.
+
+    The amplitude is randomly sampled according to the chosen distribution
+    function. `'normal'` and `'uniform'` yield continuous amplitude values and
+    are suitable to encode random amplitude and phase of reflections.
+    In contrast, `'binary'` yields only -1 and 1, and hence is only suitable
+    to encode random phase.
+
+    The final reflection sequence is generated by mapping the arrival times to
+    uniform time samples. Duplicate time samples are removed, which results in
+    a maximum of one reflection per time sample.
+
+    Parameters
+    ----------
+    arrivals : numpy.ndarray
+        Array of arrival times in seconds.
+    n_samples : int
+        Number of samples in the output sequence.
+    sampling_rate : float
+        Sampling rate in Hz.
+    distribution : Literal['uniform', 'binary', 'normal'], optional
+        Distribution of the reflection amplitudes. Default is 'normal'.
+    seed : int, None, optional
+        Seed for the random number generator. If None, a random seed is used.
+
+    Returns
+    -------
+    pyfar.Signal
+        Reflection sequence with a maximum of one reflection per time sample.
+        The sequence has a length of ``n_samples`` and a sampling rate of
+        ``sampling_rate``.
+
+    Examples
+    --------
+    Create a reflection sequence from a set of arrival times
+    and plot the result.
+
+    .. plot::
+
+        >>> import pyrato
+        >>> import numpy as np
+        >>> import pyfar as pf
+        ...
+        >>> times_of_arrival = np.asarray([.1, .3, .35, .41])
+        >>> sequence = pyrato.parametric.random_reflection_sequence(
+        ...     times_of_arrival, n_samples=50, sampling_rate=100,
+        ...     distribution='normal', seed=10)
+        >>> pf.plot.time(sequence, marker='o', linewidth=0.5)
+
+    """
+
+    rng = np.random.default_rng(seed=seed)
+
+    sample_indices = np.round(arrivals * sampling_rate).astype(int)
+    sample_indices = sample_indices[
+        (sample_indices >= 0) & (sample_indices < n_samples)]
+
+    if distribution == "normal":
+        amplitude = rng.normal(
+            loc=0,
+            scale=1,
+            size=len(sample_indices),
+        )
+    elif distribution == 'uniform':
+        # min and max values are chosen to yield unit variance
+        amplitude = rng.uniform(
+            low=-np.sqrt(3),
+            high=np.sqrt(3),
+            size=len(sample_indices),
+        )
+    elif distribution == 'binary':
+        amplitude = rng.choice(
+            [-1, 1],
+            p=[0.5, 0.5],
+            size=len(sample_indices),
+        )
+    else:
+        raise ValueError(
+            "Unknown distribution type. "
+            "Choose from 'uniform', 'binary', or 'normal'.")
+
+    sequence = np.zeros(n_samples)
+    unique_samples, unique_idx = np.unique(sample_indices, return_index=True)
+    sequence[unique_samples] = amplitude[unique_idx]
+
+    return pf.Signal(sequence, sampling_rate)
diff --git a/tests/test_parametric.py b/tests/test_parametric.py
index d24601b..49cecdd 100644
--- a/tests/test_parametric.py
+++ b/tests/test_parametric.py
@@ -7,6 +7,8 @@
 from pyrato.parametric import mean_free_path
 import pyrato as ra
 import pyrato
+import pyfar as pf
+from scipy import stats
 
 
 @pytest.mark.parametrize(("volume","reverberation_time","expected_critical_distance"),
@@ -197,3 +199,290 @@ def test_reflection_number_errors():
             times=np.linspace(0, 1, 10),
             speed_of_sound=-300,
         )
+
+
+# ======================================================================
+# _start_time_of_arrival_poisson_process
+# ======================================================================
+
+def test_start_time_of_arrival_correct_value():
+    """Return value matches the closed-form expression."""
+    volume = 100
+    speed_of_sound = 343
+    result = parametric._start_time_of_arrival_poisson_process(
+        volume, speed_of_sound)
+    expected = (
+        2 * volume * np.log(2) / (4 * np.pi * speed_of_sound**3)
+    ) ** (1 / 3)
+    npt.assert_allclose(result, expected)
+
+
+def test_start_time_of_arrival_default_speed_of_sound():
+    """Omitting speed_of_sound uses the pyfar reference value."""
+    volume = 100
+    result_default = parametric._start_time_of_arrival_poisson_process(
+        volume)
+    result_explicit = parametric._start_time_of_arrival_poisson_process(
+        volume, pf.constants.reference_speed_of_sound)
+    npt.assert_allclose(result_default, result_explicit)
+
+
+@pytest.mark.parametrize('volume', [0, -1])
+def test_start_time_of_arrival_invalid_volume(volume):
+    """Non-positive volume must raise a ValueError."""
+    with pytest.raises(ValueError, match="'volume' must be positive"):
+        parametric._start_time_of_arrival_poisson_process(volume, 343)
+
+
+@pytest.mark.parametrize('speed_of_sound', [0, -343])
+def test_start_time_of_arrival_invalid_speed_of_sound(speed_of_sound):
+    """Non-positive speed of sound must raise a ValueError."""
+    with pytest.raises(ValueError, match="speed_of_sound must be positive"):
+        parametric._start_time_of_arrival_poisson_process(100, speed_of_sound)
+
+
+# ======================================================================
+# time_of_arrival_poisson_process
+# ======================================================================
+
+def test_poisson_process_toa_kolmogorov_smirnov_statistic():
+    """
+    Test if the time of arrival intervals are drawn according to the
+    expected distribution of reflections in a room with a given volume and
+    speed of sound.
+
+    The test uses the Kolmogorov-Smirnov test to compare the empirical
+    distribution of the time of arrival intervals with the expected cumulative
+    distribution from room acoustics theory.
+    """
+    volume = 100
+    speed_of_sound = 343
+
+    times = np.linspace(0, 1, 100)
+    toa = pyrato.parametric.time_of_arrival_poisson_process(
+        volume,
+        times,
+        speed_of_sound,
+        seed=42,
+    )
+
+    def cumulative_reflections_callable(x):
+        """Normalized to fall in the range [0, 1].
+
+        All constant parameters are not relevant after normalization,
+        only the time dependency remains, which is cubic in time.
+        """
+        return (x / times[-1])**3
+
+    ks_test = stats.kstest(
+        toa,
+        cumulative_reflections_callable,
+        alternative='two-sided',
+    )
+
+    # p-value < 0.01 reject null hypothesis that the samples are drawn from
+    # the expected distribution.
+    # p-value > 0.01 fail to reject the null hypothesis that the data are not
+    # drawn from the expected distribution.
+    assert ks_test.pvalue > 0.01
+
+def test_toa_poisson_seed_reproducibility():
+    """Same seed must yield identical arrival arrays."""
+    volume = 100
+    times = np.linspace(0, 0.1, 100)
+    toa1 = parametric.time_of_arrival_poisson_process(
+        volume, times, seed=42)
+    toa2 = parametric.time_of_arrival_poisson_process(
+        volume, times, seed=42)
+    npt.assert_array_equal(toa1, toa2)
+
+
+def test_toa_poisson_arrivals_ge_t_start():
+    """All returned arrivals must be >= the Poisson-process start time."""
+    volume = 100
+    speed_of_sound = 343
+    times = np.linspace(0, 0.1, 100)
+    t_start = parametric._start_time_of_arrival_poisson_process(
+        volume, speed_of_sound)
+    toa = parametric.time_of_arrival_poisson_process(
+        volume, times, speed_of_sound, seed=0)
+    assert len(toa) > 0
+    assert np.all(toa >= t_start)
+
+
+def test_toa_poisson_arrivals_within_time_range():
+    """All returned arrivals must lie within the supplied time vector."""
+    volume = 100
+    times = np.linspace(0, 0.1, 100)
+    toa = parametric.time_of_arrival_poisson_process(volume, times, seed=0)
+    assert len(toa) > 0
+    assert np.all(toa <= times[-1])
+
+
+def test_toa_poisson_reflection_rate_limit_reduces_events():
+    """A low reflection rate limit must produce fewer events than no limit."""
+    volume = 100
+    speed_of_sound = 343
+    times = np.linspace(0, 0.1, 100)
+    seed = 0
+    toa_unlimited = parametric.time_of_arrival_poisson_process(
+        volume, times, speed_of_sound, seed=seed)
+    toa_limited = parametric.time_of_arrival_poisson_process(
+        volume, times, speed_of_sound, reflection_rate_limit=1, seed=seed)
+    assert len(toa_limited) < len(toa_unlimited)
+
+
+@pytest.mark.parametrize('volume', [0, -1])
+def test_toa_poisson_invalid_volume(volume):
+    """Non-positive volume must raise a ValueError."""
+    with pytest.raises(ValueError, match="'volume' must be positive"):
+        parametric.time_of_arrival_poisson_process(
+            volume, np.linspace(0, 1, 10))
+
+
+@pytest.mark.parametrize('speed_of_sound', [0, -343])
+def test_toa_poisson_invalid_speed_of_sound(speed_of_sound):
+    """Non-positive speed of sound must raise a ValueError."""
+    with pytest.raises(ValueError, match="speed_of_sound must be positive"):
+        parametric.time_of_arrival_poisson_process(
+            100, np.linspace(0, 1, 10), speed_of_sound=speed_of_sound)
+
+
+# ======================================================================
+# random_reflection_sequence
+# ======================================================================
+
+@pytest.mark.parametrize('distribution', ['normal', 'uniform', 'binary'])
+def test_reflection_sequence_returns_signal(distribution):
+    """Return type must be a pyfar Signal."""
+    arrivals = np.asarray([0.1, 0.3, 0.35])
+    seq = parametric.random_reflection_sequence(
+        arrivals, n_samples=50, sampling_rate=100, distribution=distribution)
+    assert isinstance(seq, pf.Signal)
+
+
+@pytest.mark.parametrize('distribution', ['normal', 'uniform', 'binary'])
+def test_reflection_sequence_length(distribution):
+    """Output signal must have exactly n_samples samples."""
+    arrivals = np.asarray([0.1, 0.3, 0.35])
+    n_samples = 50
+    seq = parametric.random_reflection_sequence(
+        arrivals, n_samples=n_samples, sampling_rate=100,
+        distribution=distribution)
+    assert seq.n_samples == n_samples
+
+
+@pytest.mark.parametrize('distribution', ['normal', 'uniform', 'binary'])
+def test_reflection_sequence_sampling_rate(distribution):
+    """Output signal must carry the requested sampling rate."""
+    arrivals = np.asarray([0.1, 0.3])
+    sampling_rate = 44100
+    seq = parametric.random_reflection_sequence(
+        arrivals, n_samples=100, sampling_rate=sampling_rate,
+        distribution=distribution)
+    assert seq.sampling_rate == sampling_rate
+
+
+@pytest.mark.parametrize('distribution', ['normal', 'uniform', 'binary'])
+def test_reflection_sequence_seed_reproducibility(distribution):
+    """Same seed must yield an identical output signal."""
+    arrivals = np.asarray([0.1, 0.3, 0.35, 0.41])
+    seq1 = parametric.random_reflection_sequence(
+        arrivals, n_samples=50, sampling_rate=100, seed=42,
+        distribution=distribution)
+    seq2 = parametric.random_reflection_sequence(
+        arrivals, n_samples=50, sampling_rate=100, seed=42,
+        distribution=distribution)
+    npt.assert_array_equal(seq1.time, seq2.time)
+
+
+@pytest.mark.parametrize('distribution', ['normal', 'uniform', 'binary'])
+def test_reflection_sequence_arrivals_out_of_range(distribution):
+    """Arrivals whose sample index >= n_samples must be ignored."""
+    # 1.0 * 100 = 100 == n_samples, so it must be excluded
+    arrivals = np.asarray([0.1, 1.0])
+    n_samples = 100
+    seq = parametric.random_reflection_sequence(
+        arrivals, n_samples=n_samples, sampling_rate=100, seed=0,
+        distribution=distribution)
+    signal = np.squeeze(seq.time)
+    assert np.count_nonzero(signal) == 1
+    assert signal[10] != 0
+
+
+@pytest.mark.parametrize('distribution', ['normal', 'uniform', 'binary'])
+def test_reflection_sequence_unique_samples(distribution):
+    """Two arrivals mapping to the same sample yield exactly one non-zero."""
+    # 0.1 and 0.1001 both round to sample index 10 at fs=100
+    arrivals = np.asarray([0.1, 0.1001])
+    seq = parametric.random_reflection_sequence(
+        arrivals, n_samples=50, sampling_rate=100, seed=0,
+        distribution=distribution)
+    assert np.count_nonzero(np.squeeze(seq.time)) == 1
+
+
+@pytest.mark.parametrize('distribution', ['normal', 'uniform', 'binary'])
+def test_reflection_sequence_nonzero_positions(distribution):
+    """Non-zero positions must equal the rounded arrival sample indices."""
+    arrivals = np.asarray([0.1, 0.3, 0.35, 0.41])
+    n_samples = 50
+    sampling_rate = 100
+    seq = parametric.random_reflection_sequence(
+        arrivals, n_samples=n_samples, sampling_rate=sampling_rate, seed=0,
+        distribution=distribution)
+    expected_indices = np.round(arrivals * sampling_rate).astype(int)
+    nonzero_indices = np.flatnonzero(np.squeeze(seq.time))
+    npt.assert_array_equal(
+        np.sort(nonzero_indices), np.sort(expected_indices))
+
+
+def test_reflection_sequence_binary_values():
+    """Binary distribution must produce values only in {-1, 0, 1}."""
+    arrivals = np.asarray([0.1, 0.2, 0.3, 0.4, 0.5])
+    seq = parametric.random_reflection_sequence(
+        arrivals, n_samples=100, sampling_rate=100,
+        distribution='binary', seed=0)
+    assert np.all(np.isin(np.squeeze(seq.time), [-1, 0, 1]))
+
+
+def test_reflection_sequence_normal_values():
+    """Normal distribution must produce continuous (non-binary) amplitudes."""
+    arrivals = np.asarray([0.1, 0.2, 0.3, 0.4, 0.5])
+    seq = parametric.random_reflection_sequence(
+        arrivals, n_samples=100, sampling_rate=100,
+        distribution='normal', seed=0)
+    nonzero = np.squeeze(seq.time)[np.squeeze(seq.time) != 0]
+    assert not np.all(np.abs(nonzero) == 1)
+
+
+def test_reflection_sequence_uniform_values():
+    """Uniform distribution must stay within [-sqrt(3), sqrt(3)]."""
+    arrivals = np.linspace(0, 0.99, 50)
+    seq = parametric.random_reflection_sequence(
+        arrivals, n_samples=100, sampling_rate=100,
+        distribution='uniform', seed=0)
+    nonzero = np.squeeze(seq.time)[np.squeeze(seq.time) != 0]
+    assert np.all(np.abs(nonzero) <= np.sqrt(3))
+
+
+def test_reflection_sequence_invalid_distribution():
+    """An unrecognised distribution name must raise a ValueError."""
+    with pytest.raises(ValueError, match="Unknown distribution"):
+        parametric.random_reflection_sequence(
+            np.asarray([0.1]), n_samples=10, sampling_rate=100,
+            distribution='invalid')
+
+
+@pytest.mark.parametrize('distribution', ['normal', 'uniform', 'binary'])
+def test_reflection_sequence_negative_arrivals_ignored(distribution):
+    """Negative arrival times must not wrap-around array boundaries."""
+    arrivals = np.asarray([-0.1, 0.1])
+    seq = parametric.random_reflection_sequence(
+        arrivals, n_samples=50, sampling_rate=100, seed=0,
+        distribution=distribution)
+    signal = np.squeeze(seq.time)
+    # Only the arrival at 0.1 s (sample 10) must be non-zero;
+    # the negative arrival must be silently dropped, not written to sample -10.
+    assert np.count_nonzero(signal) == 1
+    assert signal[10] != 0
+    assert signal[-10] == 0  # last-10th element must be untouched