Merge pull request #62 from frazane/logs-trunc

Add log score for the truncated logistic, normal and t distributions

docs/api/logarithmic.md

@@ -20,4 +20,10 @@
 
 ::: scoringrules.logs_t
 
+::: scoringrules.logs_tlogistic
+
+::: scoringrules.logs_tnormal
+
+::: scoringrules.logs_tt
+
 ::: scoringrules.logs_uniform

scoringrules/_logs.py

@@ -30,6 +30,7 @@ def logs_binomial(
     Returns
     -------
     score:
+        The LS between tLogistic(location, scale, lower, upper) and obs.
         The LS between Binomial(n, prob) and obs.
 
     Examples
@@ -368,6 +369,135 @@ def logs_t(
     return logarithmic.t(observation, df, location, scale, backend=backend)
 
 
+def logs_tlogistic(
+    observation: "ArrayLike",
+    location: "ArrayLike",
+    scale: "ArrayLike",
+    /,
+    lower: "ArrayLike" = float("-inf"),
+    upper: "ArrayLike" = float("inf"),
+    *,
+    backend: "Backend" = None,
+) -> "ArrayLike":
+    r"""Compute the logarithmic score (LS) for the truncated logistic distribution.
+
+    This score is equivalent to the negative log likelihood of the truncated logistic distribution.
+
+    Parameters
+    ----------
+    observation: ArrayLike
+        The observed values.
+    location: ArrayLike
+        Location parameter of the forecast distribution.
+    scale: ArrayLike
+        Scale parameter of the forecast distribution.
+    lower: ArrayLike
+        Lower boundary of the truncated forecast distribution.
+    upper: ArrayLike
+        Upper boundary of the truncated forecast distribution.
+
+    Returns
+    -------
+    score:
+        The LS between tLogistic(location, scale, lower, upper) and obs.
+
+    Examples
+    --------
+    >>> import scoringrules as sr
+    >>> sr.logs_tlogistic(0.0, 0.1, 0.4, -1.0, 1.0)
+    """
+    return logarithmic.tlogistic(
+        observation, location, scale, lower, upper, backend=backend
+    )
+
+
+def logs_tnormal(
+    observation: "ArrayLike",
+    location: "ArrayLike",
+    scale: "ArrayLike",
+    /,
+    lower: "ArrayLike" = float("-inf"),
+    upper: "ArrayLike" = float("inf"),
+    *,
+    backend: "Backend" = None,
+) -> "ArrayLike":
+    r"""Compute the logarithmic score (LS) for the truncated normal distribution.
+
+    This score is equivalent to the negative log likelihood of the truncated normal distribution.
+
+    Parameters
+    ----------
+    observation: ArrayLike
+        The observed values.
+    location: ArrayLike
+        Location parameter of the forecast distribution.
+    scale: ArrayLike
+        Scale parameter of the forecast distribution.
+    lower: ArrayLike
+        Lower boundary of the truncated forecast distribution.
+    upper: ArrayLike
+        Upper boundary of the truncated forecast distribution.
+
+    Returns
+    -------
+    score:
+        The LS between tNormal(location, scale, lower, upper) and obs.
+
+    Examples
+    --------
+    >>> import scoringrules as sr
+    >>> sr.logs_tnormal(0.0, 0.1, 0.4, -1.0, 1.0)
+    """
+    return logarithmic.tnormal(
+        observation, location, scale, lower, upper, backend=backend
+    )
+
+
+def logs_tt(
+    observation: "ArrayLike",
+    df: "ArrayLike",
+    /,
+    location: "ArrayLike" = 0.0,
+    scale: "ArrayLike" = 1.0,
+    lower: "ArrayLike" = float("-inf"),
+    upper: "ArrayLike" = float("inf"),
+    *,
+    backend: "Backend" = None,
+) -> "ArrayLike":
+    r"""Compute the logarithmic score (LS) for the truncated Student's t distribution.
+
+    This score is equivalent to the negative log likelihood of the truncated t distribution.
+
+    Parameters
+    ----------
+    observation: ArrayLike
+        The observed values.
+    df: ArrayLike
+        Degrees of freedom parameter of the forecast distribution.
+    location: ArrayLike
+        Location parameter of the forecast distribution.
+    scale: ArrayLike
+        Scale parameter of the forecast distribution.
+    lower: ArrayLike
+        Lower boundary of the truncated forecast distribution.
+    upper: ArrayLike
+        Upper boundary of the truncated forecast distribution.
+
+    Returns
+    -------
+    score:
+        The LS between tt(df, location, scale, lower, upper) and obs.
+
+    Examples
+    --------
+    >>> import scoringrules as sr
+    >>> sr.logs_tt(0.0, 2.0, 0.1, 0.4, -1.0, 1.0)
+    """
+    return logarithmic.tt(
+        observation, df, location, scale, lower, upper, backend=backend
+    )
+
+
 def logs_uniform(
     observation: "ArrayLike",
     min: "ArrayLike",

scoringrules/core/logarithmic.py

@@ -7,9 +7,12 @@
     _gamma_pdf,
     _hypergeo_pdf,
     _logis_pdf,
+    _logis_cdf,
     _norm_pdf,
+    _norm_cdf,
     _pois_pdf,
     _t_pdf,
+    _t_cdf,
 )
 
 if tp.TYPE_CHECKING:
@@ -85,6 +88,30 @@ def logistic(
     return -B.log(prob)
 
 
+def tlogistic(
+    obs: "ArrayLike",
+    location: "ArrayLike",
+    scale: "ArrayLike",
+    lower: "ArrayLike",
+    upper: "ArrayLike",
+    backend: "Backend" = None,
+) -> "Array":
+    """Compute the logarithmic score for the truncated logistic distribution."""
+    B = backends.active if backend is None else backends[backend]
+    obs, mu, sigma, lower, upper = map(B.asarray, (obs, location, scale, lower, upper))
+    ω = (obs - mu) / sigma
+    u = (upper - mu) / sigma
+    l = (lower - mu) / sigma
+    F_u = _logis_cdf(u, backend=backend)
+    F_l = _logis_cdf(l, backend=backend)
+    denom = F_u - F_l
+
+    ind_out = (ω < l) | (ω > u)
+    prob = _logis_pdf(ω) / sigma
+    s = B.where(ind_out, float("inf"), -B.log(prob / denom))
+    return s
+
+
 def loglogistic(
     obs: "ArrayLike",
     mulog: "ArrayLike",
@@ -163,6 +190,57 @@ def t(
     return -B.log(prob)
 
 
+def tnormal(
+    obs: "ArrayLike",
+    location: "ArrayLike",
+    scale: "ArrayLike",
+    lower: "ArrayLike",
+    upper: "ArrayLike",
+    backend: "Backend" = None,
+) -> "Array":
+    """Compute the logarithmic score for the truncated normal distribution."""
+    B = backends.active if backend is None else backends[backend]
+    obs, mu, sigma, lower, upper = map(B.asarray, (obs, location, scale, lower, upper))
+    ω = (obs - mu) / sigma
+    u = (upper - mu) / sigma
+    l = (lower - mu) / sigma
+    F_u = _norm_cdf(u, backend=backend)
+    F_l = _norm_cdf(l, backend=backend)
+    denom = F_u - F_l
+
+    ind_out = (ω < l) | (ω > u)
+    prob = _norm_pdf(ω) / sigma
+    s = B.where(ind_out, float("inf"), -B.log(prob / denom))
+    return s
+
+
+def tt(
+    obs: "ArrayLike",
+    df: "ArrayLike",
+    location: "ArrayLike",
+    scale: "ArrayLike",
+    lower: "ArrayLike",
+    upper: "ArrayLike",
+    backend: "Backend" = None,
+) -> "Array":
+    """Compute the logarithmic score for the truncated t distribution."""
+    B = backends.active if backend is None else backends[backend]
+    obs, df, mu, sigma, lower, upper = map(
+        B.asarray, (obs, df, location, scale, lower, upper)
+    )
+    ω = (obs - mu) / sigma
+    u = (upper - mu) / sigma
+    l = (lower - mu) / sigma
+    F_u = _t_cdf(u, df, backend=backend)
+    F_l = _t_cdf(l, df, backend=backend)
+    denom = F_u - F_l
+
+    ind_out = (ω < l) | (ω > u)
+    prob = _t_pdf(ω, df) / sigma
+    s = B.where(ind_out, float("inf"), -B.log(prob / denom))
+    return s
+
+
 def uniform(
     obs: "ArrayLike",
     min: "ArrayLike",

scoringrules/core/stats.py

@@ -24,6 +24,12 @@ def _logis_pdf(x: "ArrayLike", backend: "Backend" = None) -> "Array":
     return B.exp(-x) / (1 + B.exp(-x)) ** 2
 
 
+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 B.exp(-x) / (1 + B.exp(-x)) ** 2
+
+
 def _logis_cdf(x: "ArrayLike", backend: "Backend" = None) -> "Array":
     """Cumulative distribution function for the standard logistic distribution."""
     B = backends.active if backend is None else backends[backend]

tests/test_crps.py

@@ -158,11 +158,6 @@ def test_2pexponential(backend):
     expected = 6.018359
     assert np.isclose(res, expected)
 
-    obs, scale1, scale2, location = 10.5, 4.1, 0.8, 5.0
-    res0 = _crps.crps_2pexponential(obs, scale1, scale2, location, backend=backend)
-    res = _crps.crps_2pexponential(obs, scale1, scale2, location, backend=backend)
-    assert np.isclose(res, res0)
-
 
 @pytest.mark.parametrize("backend", BACKENDS)
 def test_gamma(backend):
@@ -542,24 +537,6 @@ 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

tests/test_logs.py

@@ -121,6 +121,58 @@ def test_lognormal(backend):
     assert np.isclose(res, expected)
 
 
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_tlogis(backend):
+    obs, location, scale, lower, upper = 4.9, 3.5, 2.3, 0.0, 20.0
+    res = _logs.logs_tlogistic(obs, location, scale, lower, upper, backend=backend)
+    expected = 2.11202
+    assert np.isclose(res, expected)
+
+    # aligns with logs_logistic
+    # res0 = _logs.logs_logistic(obs, location, scale, backend=backend)
+    # res = _logs.logs_tlogistic(obs, location, scale, backend=backend)
+    # assert np.isclose(res, res0)
+
+
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_tnormal(backend):
+    obs, location, scale, lower, upper = 4.2, 2.9, 2.2, 1.5, 17.3
+    res = _logs.logs_tnormal(obs, location, scale, lower, upper, backend=backend)
+    expected = 1.577806
+    assert np.isclose(res, expected)
+
+    obs, location, scale, lower, upper = -1.0, 2.9, 2.2, 1.5, 17.3
+    res = _logs.logs_tnormal(obs, location, scale, lower, upper, backend=backend)
+    expected = float("inf")
+    assert np.isclose(res, expected)
+
+    # aligns with logs_normal
+    res0 = _logs.logs_normal(obs, location, scale, backend=backend)
+    res = _logs.logs_tnormal(obs, location, scale, backend=backend)
+    assert np.isclose(res, res0)
+
+
+@pytest.mark.parametrize("backend", BACKENDS)
+def test_tt(backend):
+    if backend in ["jax", "torch", "tensorflow"]:
+        pytest.skip("Not implemented in jax, torch or tensorflow backends")
+
+    obs, df, location, scale, lower, upper = 1.9, 2.9, 3.1, 4.2, 1.5, 17.3
+    res = _logs.logs_tt(obs, df, location, scale, lower, upper, backend=backend)
+    expected = 2.002856
+    assert np.isclose(res, expected)
+
+    obs, df, location, scale, lower, upper = -1.0, 2.9, 3.1, 4.2, 1.5, 17.3
+    res = _logs.logs_tt(obs, df, location, scale, lower, upper, backend=backend)
+    expected = float("inf")
+    assert np.isclose(res, expected)
+
+    # aligns with logs_t
+    # res0 = _logs.logs_t(obs, df, location, scale, backend=backend)
+    # res = _logs.logs_tt(obs, df, location, scale, backend=backend)
+    # assert np.isclose(res, res0)
+
+
 @pytest.mark.parametrize("backend", BACKENDS)
 def test_normal(backend):
     obs, mu, sigma = 17.1, 13.8, 3.3