Add mu and sigma to LogNormal and raise AttributeError for unknown names

.pre-commit-config.yaml

@@ -42,4 +42,7 @@ repos:
         name: pylint on tests*
         language: system
         files: ^tests
-        entry: pylint --disable=invalid-name,missing-docstring,too-many-function-args,unused-variable,fixme
+        # pre-commit splits the files over parallel pylint runs, so
+        # duplicate-code only sees pairs of files that share a run. Its result
+        # depends on the CPU count and the file list, so it is disabled here.
+        entry: pylint --disable=invalid-name,missing-docstring,too-many-function-args,unused-variable,fixme,duplicate-code

RELEASE_NOTES.md

@@ -17,13 +17,15 @@
 * Set scikit-learn classifier metadata on `NGBClassifier` so that scikit-learn correctly identifies it as a classifier (PR #413)
 * Handle array-like sample weights positionally and validate them at fit time; weights with wrong length, shape, NaN, or inf now raise `ValueError` (PR #425)
 * Document saving and loading models with `skops` and test the round trip (PR #427)
+* Add `mu` and `sigma` to `LogNormal`, the mean and standard deviation of log(Y). `LogNormal.loc` and `LogNormal.scale` are deprecated and give a `FutureWarning`. `LogNormal.scale` was sigma, but `params['scale']` is exp(mu) as in `scipy.stats.lognorm`, and users mixed the two up. `params` does not change (issue #168, PR #433)
 
 ### Bug fixes
 
 * Fix the `CRPScore` metric of `Normal`, `NormalFixedMean`, and `LogNormal`, which was off by a factor of the predicted scale. Fits using `CRPScore` with these distributions now produce correct natural gradients (issue #421, PR #424)
 * Fix the scale gradient of the `Laplace` `CRPScore`, which dropped a `0.75 * scale` term (PR #419)
 * Fix censored scores and gradients for `Exponential` and `LogNormal`, which swapped censored and observed rows and had sign and extra-term errors (PR #420)
 * Fix the log-scale Fisher information of `TFixedDf` and `Cauchy`, which was off by a factor of `4 * var` (PR #422)
+* Raise `AttributeError` for an attribute that a distribution does not have. Before, the distributions built on `scipy.stats` returned `None`, so a typo such as `dist.s` failed silently. Predicted distributions can now be pickled and deep-copied (issue #168, PR #433)
 * Reset validation iteration metadata on refit so that re-calling `fit` on the same model behaves identically to fitting a fresh one (PR #426)
 * Keep `random_state` as given and seed a fresh generator in `fit`, so that `get_params` returns the original value and clones behave identically (PR #417)
 * Replace the deprecated `friedman_mse` criterion with `squared_error` for the default tree learner (issue #408, PR #409)

docs/content/1-useage.ipynb

@@ -118,7 +118,9 @@
     "| --- | --- | --- | --- |\n",
     "| `Normal` | `loc`, `scale` | `LogScore`, `CRPScore` | [`scipy.stats` normal](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.norm.html) |\n",
     "| `LogNormal` | `s`, `scale` | `LogScore`, `CRPScore` | [`scipy.stats` lognormal](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.lognorm.html) |\n",
-    "| `Exponential` | `scale` | `LogScore`, `CRPScore` | [`scipy.stats` exponential](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.expon.html) |"
+    "| `Exponential` | `scale` | `LogScore`, `CRPScore` | [`scipy.stats` exponential](https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.expon.html) |\n",
+    "\n",
+    "The parameters in this table are the keys of `params` on a predicted distribution, and they have the same meaning as in `scipy.stats`. For `LogNormal`, $\\log Y$ has a normal distribution with mean $\\mu$ and standard deviation $\\sigma$. So `params['s']` is $\\sigma$ and `params['scale']` is $e^\\mu$. The attributes `mu` and `sigma` of a `LogNormal` give $\\mu$ and $\\sigma$."
    ]
   },
   {

docs/content/5-dev.ipynb

@@ -75,11 +75,6 @@
     "    def sample(self, m):\n",
     "        return np.array([self.dist.rvs() for i in range(m)])\n",
     "    \n",
-    "    def __getattr__(self, name): # gives us access to Laplace.mean() required for RegressionDist.predict()\n",
-    "        if name in dir(self.dist):\n",
-    "            return getattr(self.dist, name)\n",
-    "        return None\n",
-    "    \n",
     "    @property\n",
     "    def params(self):\n",
     "        return {'loc':self.loc, 'scale':self.scale}"
@@ -91,7 +86,7 @@
    "source": [
     "The `fit()` method is a class method that takes a vector of observations and fits a marginal distribution. Meanwhile, `sample()` should return a $m$ samples from $P(Y|X=x)$, each of which is a vector of `len(Y)`.\n",
     "\n",
-    "Here we're taking advantage of the fact that `scipy.stats` already has the Laplace distribution implemented so we can steal its `fit()` method and put a thin wrapper around `rvs()` to get samples. We also use `__getattr__()` on the internal `scipy.stats` object to get access to its `mean()` method.\n",
+    "Here we're taking advantage of the fact that `scipy.stats` already has the Laplace distribution implemented so we can steal its `fit()` method and put a thin wrapper around `rvs()` to get samples. NGBoost gives unknown attribute names to the `scipy.stats` object in `self.dist`, so `Laplace` gets the `mean()` method of that object. Regression distributions need `mean()` for point prediction.\n",
     "\n",
     "Lastly, we write a convenience method `params()` that, when called, returns the distributional parameters as the user expects to see them, i.e. $(\\mu, b)$, not $(\\mu, \\log b)$.\n",
     "\n",
@@ -206,11 +201,6 @@
     "    def sample(self, m):\n",
     "        return np.array([self.dist.rvs() for i in range(m)])\n",
     "    \n",
-    "    def __getattr__(self, name): # gives us access to Laplace.mean() required for RegressionDist.predict()\n",
-    "        if name in dir(self.dist):\n",
-    "            return getattr(self.dist, name)\n",
-    "        return None\n",
-    "    \n",
     "    @property\n",
     "    def params(self):\n",
     "        return {'loc':self.loc, 'scale':self.scale}"

ngboost/distns/distn.py

@@ -27,6 +27,16 @@ def __getitem__(self, key):
     def __len__(self):
         return self._params.shape[1]
 
+    def __getattr__(self, name):
+        # gives access to the methods of the scipy.stats object, e.g. mean() and cdf()
+        # read dist from __dict__ so that unpickling and deepcopy do not recurse
+        dist = self.__dict__.get("dist")
+        if dist is not None and name in dir(dist):
+            return getattr(dist, name)
+        raise AttributeError(
+            f"'{type(self).__name__}' object has no attribute '{name}'"
+        )
+
     @classmethod
     def implementation(cls, Score, scores=None):
         """

ngboost/distns/exponential.py

@@ -68,11 +68,6 @@ def __init__(self, params):  # pylint: disable=super-init-not-called
         self.scale = np.exp(params[0])
         self.dist = dist(scale=self.scale)
 
-    def __getattr__(self, name):
-        if name in dir(self.dist):
-            return getattr(self.dist, name)
-        return None
-
     # should implement a `sample()` method
 
     @property

ngboost/distns/gamma.py

@@ -51,11 +51,6 @@ def fit(Y):
     def sample(self, m):
         return np.array([self.rvs() for _ in range(m)])
 
-    def __getattr__(self, name):
-        if name in dir(self.dist):
-            return getattr(self.dist, name)
-        return None
-
     @property
     def params(self):
         return {"alpha": self.alpha, "beta": self.beta}

ngboost/distns/halfnormal.py

@@ -46,11 +46,6 @@ def fit(Y):
     def sample(self, m):
         return np.array([self.dist.rvs() for i in range(m)])
 
-    def __getattr__(self, name):
-        if name in dir(self.dist):
-            return getattr(self.dist, name)
-        return None
-
     @property
     def params(self):
         return {"loc": np.zeros(shape=self.scale.shape), "scale": self.scale}

ngboost/distns/laplace.py

@@ -71,11 +71,6 @@ def fit(Y):
     def sample(self, m):
         return np.array([self.dist.rvs() for i in range(m)])
 
-    def __getattr__(self, name):
-        if name in dir(self.dist):
-            return getattr(self.dist, name)
-        return None
-
     @property
     def params(self):
         return {"loc": self.loc, "scale": self.scale}

ngboost/distns/lognormal.py

@@ -1,5 +1,7 @@
 """The NGBoost LogNormal distribution and scores"""
 
+from warnings import warn
+
 import numpy as np
 import scipy as sp
 from scipy.stats import lognorm as dist
@@ -20,23 +22,23 @@ def d_score(self, Y):
         E = Y["Event"][:, np.newaxis]
         T = Y["Time"]
         lT = np.log(T)
-        Z = (lT - self.loc) / self.scale
+        Z = (lT - self.mu) / self.sigma
 
-        D_uncens = np.zeros((self.loc.shape[0], 2))
-        D_uncens[:, 0] = (self.loc - lT) / (self.scale**2)
-        D_uncens[:, 1] = 1 - ((self.loc - lT) ** 2) / (self.scale**2)
+        D_uncens = np.zeros((self.mu.shape[0], 2))
+        D_uncens[:, 0] = (self.mu - lT) / (self.sigma**2)
+        D_uncens[:, 1] = 1 - ((self.mu - lT) ** 2) / (self.sigma**2)
 
-        D_cens = np.zeros((self.loc.shape[0], 2))
-        D_cens[:, 0] = -sp.stats.norm.pdf(lT, loc=self.loc, scale=self.scale) / (
+        D_cens = np.zeros((self.mu.shape[0], 2))
+        D_cens[:, 0] = -sp.stats.norm.pdf(lT, loc=self.mu, scale=self.sigma) / (
             1 - self.dist.cdf(T) + self.eps
         )
         D_cens[:, 1] = -Z * sp.stats.norm.pdf(Z) / (1 - self.dist.cdf(T) + self.eps)
 
         return (1 - E) * D_cens + E * D_uncens
 
     def metric(self):
-        FI = np.zeros((self.loc.shape[0], 2, 2))
-        FI[:, 0, 0] = 1 / (self.scale**2) + self.eps
+        FI = np.zeros((self.mu.shape[0], 2, 2))
+        FI[:, 0, 0] = 1 / (self.sigma**2) + self.eps
         FI[:, 1, 1] = 2
         return FI
 
@@ -46,14 +48,14 @@ def score(self, Y):
         E = Y["Event"]
         T = Y["Time"]
         lT = np.log(T)
-        Z = (lT - self.loc) / self.scale
+        Z = (lT - self.mu) / self.sigma
 
-        crps_uncens = self.scale * (
+        crps_uncens = self.sigma * (
             Z * (2 * sp.stats.norm.cdf(Z) - 1)
             + 2 * sp.stats.norm.pdf(Z)
             - 1 / np.sqrt(np.pi)
         )
-        crps_cens = self.scale * (
+        crps_cens = self.sigma * (
             Z * sp.stats.norm.cdf(Z) ** 2
             + 2 * sp.stats.norm.cdf(Z) * sp.stats.norm.pdf(Z)
             - sp.stats.norm.cdf(np.sqrt(2) * Z) / np.sqrt(np.pi)
@@ -64,29 +66,31 @@ def d_score(self, Y):
         E = Y["Event"]
         T = Y["Time"]
         lT = np.log(T)
-        Z = (lT - self.loc) / self.scale
+        Z = (lT - self.mu) / self.sigma
 
-        D = np.zeros((self.loc.shape[0], 2))
+        D = np.zeros((self.mu.shape[0], 2))
         # the censored CRPS is scale * integral of Phi(z)**2 up to Z
         D[:, 0] = E * -(2 * sp.stats.norm.cdf(Z) - 1) + (1 - E) * -(
             sp.stats.norm.cdf(Z) ** 2
         )
-        D[:, 1] = self.score(Y) + (lT - self.loc) * D[:, 0]
+        D[:, 1] = self.score(Y) + (lT - self.mu) * D[:, 0]
         return D
 
     def metric(self):
-        I = np.zeros((self.loc.shape[0], 2, 2))
+        I = np.zeros((self.mu.shape[0], 2, 2))
         I[:, 0, 0] = 2
-        I[:, 1, 1] = self.scale**2
-        I /= 2 * np.sqrt(np.pi) * self.scale[:, np.newaxis, np.newaxis]
+        I[:, 1, 1] = self.sigma**2
+        I /= 2 * np.sqrt(np.pi) * self.sigma[:, np.newaxis, np.newaxis]
         return I
 
 
 class LogNormal(RegressionDistn):
     """
     Implements the log-normal distribution for NGBoost.
 
-    The normal distribution has two parameters, s and scale (see scipy.stats.lognorm)
+    If Y is log-normal, log(Y) is normal with mean mu and standard deviation sigma.
+    The params property gives the parameters with the names of scipy.stats.lognorm:
+    s = sigma and scale = exp(mu). The attributes mu and sigma give mu and sigma.
     This distribution has both LogScore and CRPScore implemented
     for it and both work for right-censored data.
     """
@@ -97,21 +101,38 @@ class LogNormal(RegressionDistn):
     # pylint: disable=super-init-not-called
     def __init__(self, params):
         self._params = params
-        self.loc = params[0]
-        self.scale = np.exp(params[1])
-        self.dist = dist(s=self.scale, scale=np.exp(self.loc))
+        self.mu = params[0]
+        self.sigma = np.exp(params[1])
+        self.dist = dist(s=self.sigma, scale=np.exp(self.mu))
         self.eps = 1e-5
 
-    def __getattr__(self, name):
-        if name in dir(self.dist):
-            return getattr(self.dist, name)
-        return None
+    @property
+    def loc(self):
+        warn(
+            "LogNormal.loc is deprecated and will be removed in a future release. "
+            "It is mu, the mean of log(Y). Use LogNormal.mu instead.",
+            FutureWarning,
+            stacklevel=2,
+        )
+        return self.mu
+
+    @property
+    def scale(self):
+        warn(
+            "LogNormal.scale is deprecated and will be removed in a future release. "
+            "It is sigma, the standard deviation of log(Y), and it is not "
+            "params['scale'], which is exp(mu) as in scipy.stats.lognorm. "
+            "Use LogNormal.sigma or params['s'] instead.",
+            FutureWarning,
+            stacklevel=2,
+        )
+        return self.sigma
 
     # should implmenent a `sample()` method
 
     @property
     def params(self):
-        return {"s": self.scale, "scale": np.exp(self.loc)}
+        return {"s": self.sigma, "scale": np.exp(self.mu)}
 
     def fit(Y):
         m, s = sp.stats.norm.fit(np.log(Y))

ngboost/distns/normal.py

@@ -80,13 +80,6 @@ def fit(Y):
     def sample(self, m):
         return np.array([self.rvs() for i in range(m)])
 
-    def __getattr__(
-        self, name
-    ):  # gives us Normal.mean() required for RegressionDist.predict()
-        if name in dir(self.dist):
-            return getattr(self.dist, name)
-        return None
-
     @property
     def params(self):
         return {"loc": self.loc, "scale": self.scale}

ngboost/distns/poisson.py

@@ -65,13 +65,6 @@ def fit(Y):
     def sample(self, m):
         return np.array([self.dist.rvs() for i in range(m)])
 
-    def __getattr__(
-        self, name
-    ):  # gives us access to Poisson.mean() required for RegressionDist.predict()
-        if name in dir(self.dist):
-            return getattr(self.dist, name)
-        return None
-
     @property
     def params(self):
         return {"mu": self.mu}

ngboost/distns/t.py

@@ -72,9 +72,6 @@ def fit(Y):
     def sample(self, m):
         return np.array([self.rvs() for i in range(m)])
 
-    def __getattr__(self, name):
-        return getattr(self.dist, name, None)
-
     @property
     def params(self):
         return {"loc": self.loc, "scale": self.scale}
@@ -136,9 +133,6 @@ def fit(Y):
     def sample(self, m):
         return np.array([self.rvs() for i in range(m)])
 
-    def __getattr__(self, name):
-        return getattr(self.dist, name, None)
-
     @property
     def params(self):
         return {"loc": self.loc, "scale": self.scale}
@@ -194,9 +188,6 @@ def fit(Y):
     def sample(self, m):
         return np.array([self.rvs() for i in range(m)])
 
-    def __getattr__(self, name):
-        return getattr(self.dist, name, None)
-
     @property
     def params(self):
         return {"loc": self.loc, "scale": self.scale}

ngboost/distns/weibull.py

@@ -56,11 +56,6 @@ def fit(Y):
     def sample(self, m):
         return np.array([self.dist.rvs() for i in range(m)])
 
-    def __getattr__(self, name):
-        if name in dir(self.dist):
-            return getattr(self.dist, name)
-        return None
-
     @property
     def params(self):
         return {"shape": self.shape, "scale": self.scale}