From f2d4a08d0322b20953f6feb048eecb5607120085 Mon Sep 17 00:00:00 2001
From: Sam Allen <34094291+sallen12@users.noreply.github.com>
Date: Tue, 10 Sep 2024 19:50:54 +0200
Subject: [PATCH 1/2] add crps for twopiece normal distribution (#56)

* add crps for twopiece normal distribution

* fix lint issues in crps of two piece normal

* fix array bug in crps for two piece normal distribution

* fix array bug in crps for two piece normal distribution
---
 docs/api/crps.md         |  2 ++
 scoringrules/__init__.py |  2 ++
 scoringrules/_crps.py    | 67 ++++++++++++++++++++++++++++++++++++++--
 tests/test_crps.py       | 18 +++++++++++
 4 files changed, 87 insertions(+), 2 deletions(-)

diff --git a/docs/api/crps.md b/docs/api/crps.md
index 1c6a18b..ec1a589 100644
--- a/docs/api/crps.md
+++ b/docs/api/crps.md
@@ -63,6 +63,8 @@ also be viewed as the Brier score integrated over all real-valued thresholds.
 
 ::: scoringrules.crps_normal
 
+::: scoringrules.crps_2pnormal
+
 ::: scoringrules.crps_poisson
 
 ::: scoringrules.crps_t
diff --git a/scoringrules/__init__.py b/scoringrules/__init__.py
index 92b1431..709bfea 100644
--- a/scoringrules/__init__.py
+++ b/scoringrules/__init__.py
@@ -28,6 +28,7 @@
     crps_lognormal,
     crps_negbinom,
     crps_normal,
+    crps_2pnormal,
     crps_poisson,
     crps_t,
     crps_uniform,
@@ -85,6 +86,7 @@
     "crps_lognormal",
     "crps_negbinom",
     "crps_normal",
+    "crps_2pnormal",
     "crps_poisson",
     "crps_quantile",
     "crps_t",
diff --git a/scoringrules/_crps.py b/scoringrules/_crps.py
index a7addc2..8be0457 100644
--- a/scoringrules/_crps.py
+++ b/scoringrules/_crps.py
@@ -1621,12 +1621,75 @@ def crps_normal(
 
     Examples
     --------
-    >>> from scoringrules import crps
-    >>> crps.normal(0.1, 0.4, 0.0)
+    >>> import scoringrules as sr
+    >>> sr.crps_normal(0.0, 0.1, 0.4)
     """
     return crps.normal(observation, mu, sigma, backend=backend)
 
 
+def crps_2pnormal(
+    observation: "ArrayLike",
+    scale1: "ArrayLike",
+    scale2: "ArrayLike",
+    location: "ArrayLike",
+    /,
+    *,
+    backend: "Backend" = None,
+) -> "ArrayLike":
+    r"""Compute the closed form of the CRPS for the two-piece normal distribution.
+
+    It is based on the following relationship given in
+    [Jordan et al. (2019)](https://www.jstatsoft.org/article/view/v090i12):
+
+    $$ \mathrm{CRPS}(F_{\sigma_{1}, \sigma_{2}, \mu}, y) =
+    \sigma_{1} \mathrm{CRPS} \left( F_{-\infty,0}^{0, \sigma_{2}/(\sigma_{1} + \sigma_{2})}, \frac{\min(0, y - \mu)}{\sigma_{1}} \right) +
+    \sigma_{2} \mathrm{CRPS} \left( F_{0, \sigma_{1}/(\sigma_{1} + \sigma_{2})}^{\infty, 0}, \frac{\min(0, y - \mu)}{\sigma_{2}} \right), $$
+
+    where $F_{\sigma_{1}, \sigma_{2}, \mu}$ is the two-piece normal distribution with
+    scale1 and scale2 parameters $\sigma_{1}, \sigma_{2} > 0$ and location parameter $\mu$,
+    and $F_{l, L}^{u, U}$ is the CDF of the generalised truncated and censored normal distribution.
+
+    Parameters
+    ----------
+    observations: ArrayLike
+        The observed values.
+    scale1: ArrayLike
+        Scale parameter of the lower half of the forecast two-piece normal distribution.
+    scale2: ArrayLike
+        Scale parameter of the upper half of the forecast two-piece normal distribution.
+    mu: ArrayLike
+        Location parameter of the forecast two-piece normal distribution.
+
+    Returns
+    -------
+    crps: array_like
+        The CRPS between 2pNormal(scale1, scale2, mu) and obs.
+
+    Examples
+    --------
+    >>> import scoringrules as sr
+    >>> sr.crps_2pnormal(0.0, 0.4, 2.0, 0.1)
+    """
+    B = backends.active if backend is None else backends[backend]
+    lower = float("-inf")
+    upper = 0.0
+    lmass = 0.0
+    umass = scale2 / (scale1 + scale2)
+    z = B.minimum(B.asarray(0.0), B.asarray(observation - location)) / scale1
+    s1 = scale1 * crps.gtcnormal(
+        z, 0.0, 1.0, lower, upper, lmass, umass, backend=backend
+    )
+    lower = 0.0
+    upper = float("inf")
+    lmass = scale1 / (scale1 + scale2)
+    umass = 0.0
+    z = B.maximum(B.asarray(0.0), B.asarray(observation - location)) / scale2
+    s2 = scale2 * crps.gtcnormal(
+        z, 0.0, 1.0, lower, upper, lmass, umass, backend=backend
+    )
+    return s1 + s2
+
+
 def crps_poisson(
     observation: "ArrayLike",
     mean: "ArrayLike",
diff --git a/tests/test_crps.py b/tests/test_crps.py
index be5cd9d..cf354f8 100644
--- a/tests/test_crps.py
+++ b/tests/test_crps.py
@@ -542,6 +542,24 @@ def test_normal(backend):
     assert not np.any(res - 0.0 > 0.0001)
 
 
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_2pnormal(backend):
+    obs, scale1, scale2, location = 29.1, 4.6, 1.3, 27.9
+    expected = 2.189609
+    res = _crps.crps_2pnormal(obs, scale1, scale2, location, backend=backend)
+    assert np.isclose(res, expected)
+
+    obs, scale1, scale2, location = -2.2, 1.6, 3.3, -1.9
+    expected = 0.8979951
+    res = _crps.crps_2pnormal(obs, scale1, scale2, location, backend=backend)
+    assert np.isclose(res, expected)
+
+    obs, scale, location = 1.5, 4.5, 5.4
+    res0 = _crps.crps_normal(obs, location, scale, backend=backend)
+    res = _crps.crps_2pnormal(obs, scale, scale, location, backend=backend)
+    assert np.isclose(res, res0)
+
+
 @pytest.mark.parametrize("backend", BACKENDS)
 def test_poisson(backend):
     obs, mean = 1.0, 3.0

From cfe1ec947a8e2e69ca9329afa2379ad571eb7065 Mon Sep 17 00:00:00 2001
From: Sam Allen <34094291+sallen12@users.noreply.github.com>
Date: Tue, 10 Sep 2024 19:53:45 +0200
Subject: [PATCH 2/2] fix lint issues in log scores (#61)

---
 docs/api/logarithmic.md          |  20 ++
 scoringrules/_crps.py            |   2 +-
 scoringrules/_logs.py            | 364 ++++++++++++++++++++++++++++++-
 scoringrules/core/logarithmic.py | 157 ++++++++++++-
 scoringrules/core/stats.py       |  43 +++-
 tests/test_logs.py               | 180 ++++++++++++++-
 6 files changed, 744 insertions(+), 22 deletions(-)

diff --git a/docs/api/logarithmic.md b/docs/api/logarithmic.md
index e803747..4a4e737 100644
--- a/docs/api/logarithmic.md
+++ b/docs/api/logarithmic.md
@@ -1,3 +1,23 @@
 # Logarithmic Score
 
+::: scoringrules.logs_binomial
+
+::: scoringrules.logs_exponential
+
+::: scoringrules.logs_gamma
+
+::: scoringrules.logs_hypergeometric
+
+::: scoringrules.logs_logistic
+
+::: scoringrules.logs_loglogistic
+
+::: scoringrules.logs_lognormal
+
 ::: scoringrules.logs_normal
+
+::: scoringrules.logs_poisson
+
+::: scoringrules.logs_t
+
+::: scoringrules.logs_uniform
diff --git a/scoringrules/_crps.py b/scoringrules/_crps.py
index 8be0457..ae8979f 100644
--- a/scoringrules/_crps.py
+++ b/scoringrules/_crps.py
@@ -1712,7 +1712,7 @@ def crps_poisson(
     observation: ArrayLike
         The observed values.
     mean: ArrayLike
-        Mean parameter of the forecast exponential distribution.
+        Mean parameter of the forecast poisson distribution.
 
     Returns
     -------
diff --git a/scoringrules/_logs.py b/scoringrules/_logs.py
index 2298f57..f2e7f15 100644
--- a/scoringrules/_logs.py
+++ b/scoringrules/_logs.py
@@ -6,6 +6,261 @@
     from scoringrules.core.typing import Array, ArrayLike, Backend
 
 
+def logs_binomial(
+    observation: "ArrayLike",
+    n: "ArrayLike",
+    prob: "ArrayLike",
+    /,
+    *,
+    backend: "Backend" = None,
+) -> "ArrayLike":
+    r"""Compute the logarithmic score (LS) for the binomial distribution.
+
+    This score is equivalent to the negative log likelihood of the binomial distribution
+
+    Parameters
+    ----------
+    observation:
+        The observed values as an integer or array of integers.
+    n:
+        Size parameter of the forecast binomial distribution as an integer or array of integers.
+    prob:
+        Probability parameter of the forecast binomial distribution as a float or array of floats.
+
+    Returns
+    -------
+    score:
+        The LS between Binomial(n, prob) and obs.
+
+    Examples
+    --------
+    >>> import scoringrules as sr
+    >>> sr.logs_binomial(4, 10, 0.5)
+    """
+    return logarithmic.binomial(observation, n, prob, backend=backend)
+
+
+def logs_exponential(
+    observation: "ArrayLike",
+    rate: "ArrayLike",
+    /,
+    *,
+    backend: "Backend" = None,
+) -> "ArrayLike":
+    r"""Compute the logarithmic score (LS) for the exponential distribution.
+
+    This score is equivalent to the negative log likelihood of the exponential distribution
+
+    Parameters
+    ----------
+    observation:
+        The observed values.
+    rate:
+        Rate parameter of the forecast exponential distribution.
+
+    Returns
+    -------
+    score:
+        The LS between Exp(rate) and obs.
+
+    Examples
+    --------
+    >>> import scoringrules as sr
+    >>> sr.logs_exponential(0.8, 3.0)
+    """
+    return logarithmic.exponential(observation, rate, backend=backend)
+
+
+def logs_gamma(
+    observation: "ArrayLike",
+    shape: "ArrayLike",
+    /,
+    rate: "ArrayLike | None" = None,
+    *,
+    scale: "ArrayLike | None" = None,
+    backend: "Backend" = None,
+) -> "ArrayLike":
+    r"""Compute the logarithmic score (LS) for the gamma distribution.
+
+    This score is equivalent to the negative log likelihood of the gamma distribution
+
+    Parameters
+    ----------
+    observation:
+        The observed values.
+    shape:
+        Shape parameter of the forecast gamma distribution.
+    rate:
+        Rate parameter of the forecast gamma distribution.
+    scale:
+        Scale parameter of the forecast gamma distribution, where `scale = 1 / rate`.
+
+    Returns
+    -------
+    score:
+        The LS between obs and Gamma(shape, rate).
+
+    Examples
+    --------
+    >>> import scoringrules as sr
+    >>> sr.logs_gamma(0.2, 1.1, 0.1)
+
+    Raises
+    ------
+    ValueError
+        If both `rate` and `scale` are provided, or if neither is provided.
+    """
+    if (scale is None and rate is None) or (scale is not None and rate is not None):
+        raise ValueError(
+            "Either `rate` or `scale` must be provided, but not both or neither."
+        )
+
+    if rate is None:
+        rate = 1.0 / scale
+
+    return logarithmic.gamma(observation, shape, rate, backend=backend)
+
+
+def logs_hypergeometric(
+    observation: "ArrayLike",
+    m: "ArrayLike",
+    n: "ArrayLike",
+    k: "ArrayLike",
+    /,
+    *,
+    backend: "Backend" = None,
+) -> "ArrayLike":
+    r"""Compute the logarithmic score (LS) for the hypergeometric distribution.
+
+    This score is equivalent to the negative log likelihood of the hypergeometric distribution
+
+    Parameters
+    ----------
+    observation:
+        The observed values.
+    m:
+        Number of success states in the population.
+    n:
+        Number of failure states in the population.
+    k:
+        Number of draws, without replacement. Must be in 0, 1, ..., m + n.
+    backend:
+        The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'.
+
+    Returns
+    -------
+    score:
+        The LS between obs and Hypergeometric(m, n, k).
+
+    Examples
+    --------
+    >>> import scoringrules as sr
+    >>> sr.logs_hypergeometric(5, 7, 13, 12)
+    """
+    return logarithmic.hypergeometric(observation, m, n, k, backend=backend)
+
+
+def logs_logistic(
+    observation: "ArrayLike",
+    mu: "ArrayLike",
+    sigma: "ArrayLike",
+    /,
+    *,
+    backend: "Backend" = None,
+) -> "ArrayLike":
+    r"""Compute the logarithmic score (LS) for the logistic distribution.
+
+    This score is equivalent to the negative log likelihood of the logistic distribution
+
+    Parameters
+    ----------
+    observations: ArrayLike
+        Observed values.
+    mu: ArrayLike
+        Location parameter of the forecast logistic distribution.
+    sigma: ArrayLike
+        Scale parameter of the forecast logistic distribution.
+
+    Returns
+    -------
+    score:
+        The LS for the Logistic(mu, sigma) forecasts given the observations.
+
+    Examples
+    --------
+    >>> import scoringrules as sr
+    >>> sr.logs_logistic(0.0, 0.4, 0.1)
+    """
+    return logarithmic.logistic(observation, mu, sigma, backend=backend)
+
+
+def logs_loglogistic(
+    observation: "ArrayLike",
+    mulog: "ArrayLike",
+    sigmalog: "ArrayLike",
+    backend: "Backend" = None,
+) -> "ArrayLike":
+    r"""Compute the logarithmic score (LS) for the log-logistic distribution.
+
+    This score is equivalent to the negative log likelihood of the log-logistic distribution
+
+    Parameters
+    ----------
+    observation:
+        The observed values.
+    mulog:
+        Location parameter of the log-logistic distribution.
+    sigmalog:
+        Scale parameter of the log-logistic distribution.
+    backend:
+        The name of the backend used for computations. Defaults to 'numba' if available, else 'numpy'.
+
+
+    Returns
+    -------
+    score:
+        The LS between obs and Loglogis(mulog, sigmalog).
+
+    Examples
+    --------
+    >>> import scoringrules as sr
+    >>> sr.logs_loglogistic(3.0, 0.1, 0.9)
+    """
+    return logarithmic.loglogistic(observation, mulog, sigmalog, backend=backend)
+
+
+def logs_lognormal(
+    observation: "ArrayLike",
+    mulog: "ArrayLike",
+    sigmalog: "ArrayLike",
+    backend: "Backend" = None,
+) -> "ArrayLike":
+    r"""Compute the logarithmic score (LS) for the log-normal distribution.
+
+    This score is equivalent to the negative log likelihood of the log-normal distribution
+
+    Parameters
+    ----------
+    observation:
+        The observed values.
+    mulog:
+        Mean of the normal underlying distribution.
+    sigmalog:
+        Standard deviation of the underlying normal distribution.
+
+    Returns
+    -------
+    score:
+        The LS between Lognormal(mu, sigma) and obs.
+
+    Examples
+    --------
+    >>> import scoringrules as sr
+    >>> sr.logs_lognormal(0.0, 0.4, 0.1)
+    """
+    return logarithmic.lognormal(observation, mulog, sigmalog, backend=backend)
+
+
 def logs_normal(
     observation: "ArrayLike",
     mu: "ArrayLike",
@@ -17,7 +272,7 @@ def logs_normal(
 ) -> "Array":
     r"""Compute the logarithmic score (LS) for the normal distribution.
 
-    This score is equivalent to the negative log likelihood (if `negative = True`)
+    This score is equivalent to the (negative) log likelihood (if `negative = True`)
 
     Parameters
     ----------
@@ -32,15 +287,116 @@ def logs_normal(
 
     Returns
     -------
-    ls: array_like
+    score:
         The LS between Normal(mu, sigma) and obs.
 
     Examples
     --------
     >>> import scoringrules as sr
-    >>> sr.logs_normal(0.1, 0.4, 0.0)
-    >>> 0.033898
+    >>> sr.logs_normal(0.0, 0.4, 0.1)
     """
     return logarithmic.normal(
         observation, mu, sigma, negative=negative, backend=backend
     )
+
+
+def logs_poisson(
+    observation: "ArrayLike",
+    mean: "ArrayLike",
+    /,
+    *,
+    backend: "Backend" = None,
+) -> "ArrayLike":
+    r"""Compute the logarithmic score (LS) for the Poisson distribution.
+
+    This score is equivalent to the negative log likelihood of the Poisson distribution.
+
+    Parameters
+    ----------
+    observation: ArrayLike
+        The observed values.
+    mean: ArrayLike
+        Mean parameter of the forecast poisson distribution.
+
+    Returns
+    -------
+    score:
+        The LS between Pois(mean) and obs.
+
+    Examples
+    --------
+    >>> import scoringrules as sr
+    >>> sr.logs_poisson(1, 2)
+    """
+    return logarithmic.poisson(observation, mean, backend=backend)
+
+
+def logs_t(
+    observation: "ArrayLike",
+    df: "ArrayLike",
+    /,
+    location: "ArrayLike" = 0.0,
+    scale: "ArrayLike" = 1.0,
+    *,
+    backend: "Backend" = None,
+) -> "ArrayLike":
+    r"""Compute the logarithmic score (LS) for the Student's t distribution.
+
+    This score is equivalent to the negative log likelihood of the t distribution.
+
+    Parameters
+    ----------
+    observation: ArrayLike
+        The observed values.
+    df: ArrayLike
+        Degrees of freedom parameter of the forecast t distribution.
+    location: ArrayLike
+        Location parameter of the forecast t distribution.
+    sigma: ArrayLike
+        Scale parameter of the forecast t distribution.
+
+    Returns
+    -------
+    score:
+        The LS between t(df, location, scale) and obs.
+
+    Examples
+    --------
+    >>> import scoringrules as sr
+    >>> sr.logs_t(0.0, 0.1, 0.4, 0.1)
+    """
+    return logarithmic.t(observation, df, location, scale, backend=backend)
+
+
+def logs_uniform(
+    observation: "ArrayLike",
+    min: "ArrayLike",
+    max: "ArrayLike",
+    /,
+    *,
+    backend: "Backend" = None,
+) -> "ArrayLike":
+    r"""Compute the logarithmic score (LS) for the uniform distribution.
+
+    This score is equivalent to the negative log likelihood of the uniform distribution.
+
+    Parameters
+    ----------
+    observation: ArrayLike
+        The observed values.
+    min: ArrayLike
+        Lower bound of the forecast uniform distribution.
+    max: ArrayLike
+        Upper bound of the forecast uniform distribution.
+
+    Returns
+    -------
+    score:
+        The LS between U(min, max, lmass, umass) and obs.
+
+    Examples
+    --------
+    >>> import scoringrules as sr
+    >>> sr.logs_uniform(0.4, 0.0, 1.0)
+    """
+    return logarithmic.uniform(observation, min, max, backend=backend)
diff --git a/scoringrules/core/logarithmic.py b/scoringrules/core/logarithmic.py
index 4cd5982..93b5f24 100644
--- a/scoringrules/core/logarithmic.py
+++ b/scoringrules/core/logarithmic.py
@@ -1,12 +1,124 @@
 import typing as tp
 
 from scoringrules.backend import backends
-from scoringrules.core.stats import _norm_pdf
+from scoringrules.core.stats import (
+    _binom_pdf,
+    _exp_pdf,
+    _gamma_pdf,
+    _hypergeo_pdf,
+    _logis_pdf,
+    _norm_pdf,
+    _pois_pdf,
+    _t_pdf,
+)
 
 if tp.TYPE_CHECKING:
     from scoringrules.core.typing import Array, ArrayLike, Backend
 
 
+def binomial(
+    obs: "ArrayLike",
+    n: "ArrayLike",
+    prob: "ArrayLike",
+    backend: "Backend" = None,
+) -> "Array":
+    """Compute the logarithmic score for the binomial distribution."""
+    B = backends.active if backend is None else backends[backend]
+    obs, n, prob = map(B.asarray, (obs, n, prob))
+    zind = obs < 0.0
+    obs = B.where(zind, B.nan, obs)
+    prob = _binom_pdf(obs, n, prob, backend=backend)
+    s = B.where(zind, float("inf"), -B.log(prob))
+    return s
+
+
+def exponential(
+    obs: "ArrayLike", rate: "ArrayLike", backend: "Backend" = None
+) -> "Array":
+    """Compute the logarithmic score for the exponential distribution."""
+    B = backends.active if backend is None else backends[backend]
+    rate, obs = map(B.asarray, (rate, obs))
+    zind = obs < 0.0
+    obs = B.where(zind, B.nan, obs)
+    prob = _exp_pdf(obs, rate, backend=backend)
+    s = B.where(zind, float("inf"), -B.log(prob))
+    return s
+
+
+def gamma(
+    obs: "ArrayLike", shape: "ArrayLike", rate: "ArrayLike", backend: "Backend" = None
+) -> "Array":
+    """Compute the logarithmic score for the gamma distribution."""
+    B = backends.active if backend is None else backends[backend]
+    obs, shape, rate = map(B.asarray, (obs, shape, rate))
+    prob = _gamma_pdf(obs, shape, rate, backend=backend)
+    return -B.log(prob)
+
+
+def hypergeometric(
+    obs: "ArrayLike",
+    m: "ArrayLike",
+    n: "ArrayLike",
+    k: "ArrayLike",
+    backend: "Backend" = None,
+) -> "Array":
+    """Compute the logarithmic score for the hypergeometric distribution."""
+    B = backends.active if backend is None else backends[backend]
+    obs, m, n, k = map(B.asarray, (obs, m, n, k))
+    M = m + n
+    N = k
+    prob = _hypergeo_pdf(obs, M, m, N, backend=backend)
+    return -B.log(prob)
+
+
+def logistic(
+    obs: "ArrayLike",
+    mu: "ArrayLike",
+    sigma: "ArrayLike",
+    backend: "Backend" = None,
+) -> "Array":
+    """Compute the logarithmic score for the logistic distribution."""
+    B = backends.active if backend is None else backends[backend]
+    mu, sigma, obs = map(B.asarray, (mu, sigma, obs))
+    ω = (obs - mu) / sigma
+    prob = _logis_pdf(ω, backend=backend) / sigma
+    return -B.log(prob)
+
+
+def loglogistic(
+    obs: "ArrayLike",
+    mulog: "ArrayLike",
+    sigmalog: "ArrayLike",
+    backend: "Backend" = None,
+) -> "Array":
+    """Compute the logarithmic score for the log-logistic distribution."""
+    B = backends.active if backend is None else backends[backend]
+    mulog, sigmalog, obs = map(B.asarray, (mulog, sigmalog, obs))
+    zind = obs <= 0.0
+    obs = B.where(zind, B.nan, obs)
+    ω = (B.log(obs) - mulog) / sigmalog
+    prob = _logis_pdf(ω, backend=backend) / sigmalog
+    s = B.where(zind, float("inf"), -B.log(prob / obs))
+    return s
+
+
+def lognormal(
+    obs: "ArrayLike",
+    mulog: "ArrayLike",
+    sigmalog: "ArrayLike",
+    backend: "Backend" = None,
+) -> "Array":
+    """Compute the logarithmic score for the log-normal distribution."""
+    B = backends.active if backend is None else backends[backend]
+    mulog, sigmalog, obs = map(B.asarray, (mulog, sigmalog, obs))
+    zind = obs <= 0.0
+    obs = B.where(zind, B.nan, obs)
+    ω = (B.log(obs) - mulog) / sigmalog
+    prob = _norm_pdf(ω, backend=backend) / sigmalog
+    s = B.where(zind, float("inf"), -B.log(prob / obs))
+    return s
+
+
 def normal(
     obs: "ArrayLike",
     mu: "ArrayLike",
@@ -20,5 +132,46 @@ def normal(
 
     constant = -1.0 if negative else 1.0
     ω = (obs - mu) / sigma
-    prob = _norm_pdf(ω) / sigma
+    prob = _norm_pdf(ω, backend=backend) / sigma
     return constant * B.log(prob)
+
+
+def poisson(
+    obs: "ArrayLike",
+    mean: "ArrayLike",
+    backend: "Backend" = None,
+) -> "Array":
+    """Compute the logarithmic score for the Poisson distribution."""
+    B = backends.active if backend is None else backends[backend]
+    mean, obs = map(B.asarray, (mean, obs))
+    prob = _pois_pdf(obs, mean, backend=backend)
+    return -B.log(prob)
+
+
+def t(
+    obs: "ArrayLike",
+    df: "ArrayLike",
+    location: "ArrayLike",
+    scale: "ArrayLike",
+    backend: "Backend" = None,
+) -> "Array":
+    """Compute the logarithmic score for the t distribution."""
+    B = backends.active if backend is None else backends[backend]
+    df, mu, sigma, obs = map(B.asarray, (df, location, scale, obs))
+    ω = (obs - mu) / sigma
+    prob = _t_pdf(ω, df, backend=backend) / sigma
+    return -B.log(prob)
+
+
+def uniform(
+    obs: "ArrayLike",
+    min: "ArrayLike",
+    max: "ArrayLike",
+    backend: "Backend" = None,
+) -> "Array":
+    """Compute the logarithmic score for the uniform distribution."""
+    B = backends.active if backend is None else backends[backend]
+    min, max, obs = map(B.asarray, (min, max, obs))
+    ind_out = (obs > max) | (obs < min)
+    prob = B.where(ind_out, 0.0, 1.0 / (max - min))
+    return -B.log(prob)
diff --git a/scoringrules/core/stats.py b/scoringrules/core/stats.py
index 881fb22..1153dc6 100644
--- a/scoringrules/core/stats.py
+++ b/scoringrules/core/stats.py
@@ -6,16 +6,22 @@
     from scoringrules.core.typing import Array, ArrayLike, Backend
 
 
+def _norm_pdf(x: "ArrayLike", backend: "Backend" = None) -> "Array":
+    """Probability density function for the standard normal distribution."""
+    B = backends.active if backend is None else backends[backend]
+    return (1.0 / B.sqrt(2.0 * B.pi)) * B.exp(-(x**2) / 2)
+
+
 def _norm_cdf(x: "ArrayLike", backend: "Backend" = None) -> "Array":
     """Cumulative distribution function for the standard normal distribution."""
     B = backends.active if backend is None else backends[backend]
     return (1.0 + B.erf(x / B.sqrt(2.0))) / 2.0
 
 
-def _norm_pdf(x: "ArrayLike", backend: "Backend" = None) -> "Array":
-    """Probability density function for the standard normal distribution."""
+def _logis_pdf(x: "ArrayLike", backend: "Backend" = None) -> "Array":
+    """Probability density function for the standard logistic distribution."""
     B = backends.active if backend is None else backends[backend]
-    return (1.0 / B.sqrt(2.0 * B.pi)) * B.exp(-(x**2) / 2)
+    return B.exp(-x) / (1 + B.exp(-x)) ** 2
 
 
 def _logis_cdf(x: "ArrayLike", backend: "Backend" = None) -> "Array":
@@ -24,12 +30,27 @@ def _logis_cdf(x: "ArrayLike", backend: "Backend" = None) -> "Array":
     return 1 / (1 + B.exp(-x))
 
 
+def _exp_pdf(x: "ArrayLike", rate: "ArrayLike", backend: "Backend" = None) -> "Array":
+    """Probability density function for the exponential distribution."""
+    B = backends.active if backend is None else backends[backend]
+    return B.where(x < 0.0, 0.0, rate * B.exp(-rate * x))
+
+
 def _exp_cdf(x: "ArrayLike", rate: "ArrayLike", backend: "Backend" = None) -> "Array":
     """Cumulative distribution function for the exponential distribution."""
     B = backends.active if backend is None else backends[backend]
     return B.maximum(1 - B.exp(-rate * x), B.asarray(0.0))
 
 
+def _gamma_pdf(
+    x: "ArrayLike", shape: "ArrayLike", rate: "ArrayLike", backend: "Backend" = None
+) -> "Array":
+    """Probability density function for the gamma distribution."""
+    B = backends.active if backend is None else backends[backend]
+    prob = (rate**shape) * (x ** (shape - 1)) * (B.exp(-rate * x)) / B.gamma(shape)
+    return B.where(x <= 0.0, 0.0, prob)
+
+
 def _gamma_cdf(
     x: "ArrayLike", shape: "ArrayLike", rate: "ArrayLike", backend: "Backend" = None
 ) -> "Array":
@@ -39,14 +60,6 @@ def _gamma_cdf(
     return B.maximum(B.gammainc(shape, rate * B.maximum(x, zero)), zero)
 
 
-def _pois_cdf(x: "ArrayLike", mean: "ArrayLike", backend: "Backend" = None) -> "Array":
-    """Cumulative distribution function for the Poisson distribution."""
-    B = backends.active if backend is None else backends[backend]
-    x_plus = B.abs(x)
-    p = B.gammauinc(B.floor(x_plus + 1), mean) / B.gamma(B.floor(x_plus + 1))
-    return B.where(x < 0.0, 0.0, p)
-
-
 def _pois_pdf(x: "ArrayLike", mean: "ArrayLike", backend: "Backend" = None) -> "Array":
     """Probability mass function for the Poisson distribution."""
     B = backends.active if backend is None else backends[backend]
@@ -59,6 +72,14 @@ def _pois_pdf(x: "ArrayLike", mean: "ArrayLike", backend: "Backend" = None) -> "
     return B.where(mean < 0.0, B.nan, B.where(x < 0.0, 0.0, d))
 
 
+def _pois_cdf(x: "ArrayLike", mean: "ArrayLike", backend: "Backend" = None) -> "Array":
+    """Cumulative distribution function for the Poisson distribution."""
+    B = backends.active if backend is None else backends[backend]
+    x_plus = B.abs(x)
+    p = B.gammauinc(B.floor(x_plus + 1), mean) / B.gamma(B.floor(x_plus + 1))
+    return B.where(x < 0.0, 0.0, p)
+
+
 def _t_pdf(x: "ArrayLike", df: "ArrayLike", backend: "Backend" = None) -> "Array":
     """Probability density function for the standard Student's t distribution."""
     B = backends.active if backend is None else backends[backend]
diff --git a/tests/test_logs.py b/tests/test_logs.py
index 4b0706b..30d7d6b 100644
--- a/tests/test_logs.py
+++ b/tests/test_logs.py
@@ -4,11 +4,183 @@
 
 from .conftest import BACKENDS
 
-ENSEMBLE_SIZE = 51
-N = 100
+
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_binomial(backend):
+    # test correctness
+    res = _logs.logs_binomial(8, 10, 0.9, backend=backend)
+    expected = 1.641392
+    assert np.isclose(res, expected)
+
+    res = _logs.logs_binomial(-8, 10, 0.9, backend=backend)
+    expected = float("inf")
+    assert np.isclose(res, expected)
+
+    res = _logs.logs_binomial(18, 30, 0.5, backend=backend)
+    expected = 2.518839
+    assert np.isclose(res, expected)
+
+
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_exponential(backend):
+    obs, rate = 0.3, 0.1
+    res = _logs.logs_exponential(obs, rate, backend=backend)
+    expected = 2.332585
+    assert np.isclose(res, expected)
+
+    obs, rate = -1.3, 2.4
+    res = _logs.logs_exponential(obs, rate, backend=backend)
+    expected = float("inf")
+    assert np.isclose(res, expected)
+
+    obs, rate = 0.0, 0.9
+    res = _logs.logs_exponential(obs, rate, backend=backend)
+    expected = 0.1053605
+    assert np.isclose(res, expected)
+
+
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_gamma(backend):
+    obs, shape, rate = 0.2, 1.1, 0.7
+    expected = 0.6434138
+
+    res = _logs.logs_gamma(obs, shape, rate, backend=backend)
+    assert np.isclose(res, expected)
+
+    res = _logs.logs_gamma(obs, shape, scale=1 / rate, backend=backend)
+    assert np.isclose(res, expected)
+
+    with pytest.raises(ValueError):
+        _logs.logs_gamma(obs, shape, rate, scale=1 / rate, backend=backend)
+        return
+
+    with pytest.raises(ValueError):
+        _logs.logs_gamma(obs, shape, backend=backend)
+        return
+
+
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_hypergeometric(backend):
+    res = _logs.logs_hypergeometric(5, 7, 13, 12)
+    expected = 1.251525
+    assert np.isclose(res, expected)
+
+    res = _logs.logs_hypergeometric(5 * np.ones((2, 2)), 7, 13, 12, backend=backend)
+    assert res.shape == (2, 2)
+
+    res = _logs.logs_hypergeometric(5, 7 * np.ones((2, 2)), 13, 12, backend=backend)
+    assert res.shape == (2, 2)
+
+
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_logis(backend):
+    obs, mu, sigma = 17.1, 13.8, 3.3
+    res = _logs.logs_logistic(obs, mu, sigma, backend=backend)
+    expected = 2.820446
+    assert np.isclose(res, expected)
+
+    obs, mu, sigma = 3.1, 4.0, 0.5
+    res = _logs.logs_logistic(obs, mu, sigma, backend=backend)
+    expected = 1.412808
+    assert np.isclose(res, expected)
+
+
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_loglogistic(backend):
+    obs, mulog, sigmalog = 3.0, 0.1, 0.9
+    res = _logs.logs_loglogistic(obs, mulog, sigmalog, backend=backend)
+    expected = 2.672729
+    assert np.isclose(res, expected)
+
+    obs, mulog, sigmalog = 0.0, 0.1, 0.9
+    res = _logs.logs_loglogistic(obs, mulog, sigmalog, backend=backend)
+    expected = float("inf")
+    assert np.isclose(res, expected)
+
+    obs, mulog, sigmalog = 12.0, 12.1, 4.9
+    res = _logs.logs_loglogistic(obs, mulog, sigmalog, backend=backend)
+    expected = 6.299409
+    assert np.isclose(res, expected)
+
+
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_lognormal(backend):
+    obs, mulog, sigmalog = 3.0, 0.1, 0.9
+    res = _logs.logs_lognormal(obs, mulog, sigmalog, backend=backend)
+    expected = 2.527762
+    assert np.isclose(res, expected)
+
+    obs, mulog, sigmalog = 0.0, 0.1, 0.9
+    res = _logs.logs_lognormal(obs, mulog, sigmalog, backend=backend)
+    expected = float("inf")
+    assert np.isclose(res, expected)
+
+    obs, mulog, sigmalog = 12.0, 12.1, 4.9
+    res = _logs.logs_lognormal(obs, mulog, sigmalog, backend=backend)
+    expected = 6.91832
+    assert np.isclose(res, expected)
 
 
 @pytest.mark.parametrize("backend", BACKENDS)
 def test_normal(backend):
-    res = _logs.logs_normal(0.0, 0.1, 0.1, backend=backend)
-    assert np.isclose(res, -0.8836466, rtol=1e-5)
+    obs, mu, sigma = 17.1, 13.8, 3.3
+    res = _logs.logs_normal(obs, mu, sigma, backend=backend)
+    expected = 2.612861
+    assert np.isclose(res, expected)
+
+    obs, mu, sigma = 3.1, 4.0, 0.5
+    res = _logs.logs_normal(obs, mu, sigma, backend=backend)
+    expected = 1.845791
+    assert np.isclose(res, expected)
+
+
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_poisson(backend):
+    obs, mean = 1.0, 3.0
+    res = _logs.logs_poisson(obs, mean, backend=backend)
+    expected = 1.901388
+    assert np.isclose(res, expected)
+
+    obs, mean = 1.5, 2.3
+    res = _logs.logs_poisson(obs, mean, backend=backend)
+    expected = float("inf")
+    assert np.isclose(res, expected)
+
+    obs, mean = -1.0, 1.5
+    res = _logs.logs_poisson(obs, mean, backend=backend)
+    expected = float("inf")
+    assert np.isclose(res, expected)
+
+
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_t(backend):
+    if backend in ["jax", "torch", "tensorflow"]:
+        pytest.skip("Not implemented in jax, torch or tensorflow backends")
+
+    obs, df, mu, sigma = 11.1, 5.2, 13.8, 2.3
+    res = _logs.logs_t(obs, df, mu, sigma, backend=backend)
+    expected = 2.528398
+    assert np.isclose(res, expected)
+
+    obs, df = 0.7, 4.0
+    res = _logs.logs_t(obs, df, backend=backend)
+    expected = 1.269725
+    assert np.isclose(res, expected)
+
+
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_uniform(backend):
+    obs, min, max = 0.3, -1.0, 2.1
+    res = _logs.logs_uniform(obs, min, max, backend=backend)
+    expected = 1.131402
+    assert np.isclose(res, expected)
+
+    obs, min, max = -17.9, -15.2, -8.7
+    res = _logs.logs_uniform(obs, min, max, backend=backend)
+    expected = float("inf")
+    assert np.isclose(res, expected)
+
+    obs, min, max = 0.1, 0.1, 3.1
+    res = _logs.logs_uniform(obs, min, max, backend=backend)
+    expected = 1.098612
+    assert np.isclose(res, expected)