Make fit an instance method and deprecate fit_instance (#53)

* Make fit an instance method and deprecate fit_instance

* Update documentation

Author: jonathanberthias

docs/source/modules/distributions.rst

@@ -23,7 +23,6 @@ some parameters have been adjusted (e.g., the GEV distribution) to align with st
 
    ~Distribution.cdf
    ~Distribution.fit
-   ~Distribution.fit_instance
    ~Distribution.inverse_cdf
    ~Distribution.isf
    ~Distribution.logcdf

docs/source/user_guide/fitting.rst

@@ -14,16 +14,15 @@ Let's generate a larger sample from our previous object:
 
 We can fit a ``Normal`` distribution to this data, which will return another ``Normal`` object:
 
->>> Normal.fit(data)
+>>> Normal().fit(data)
 Normal(loc=1.0250822420920338, scale=1.9376400770300832)
 
 As you can see, the values are slightly different from the moments in the data.
 This is due to the fact that the ``fit`` method returns the Maximum Likelihood Estimator (MLE)
-for the data, and is thus the result of an optimisation (using **scipy.optimize**). Custom optimizer and arguments passed
-to ``scipy.optimize.minimize`` can be passed as ``kwargs`` to the ``fit`` method of any distribution.
+for the data, and is thus the result of an optimisation (using **scipy.optimize**).
 
-The syntax ``distribution.fit(data, loc=0)`` can be used to fit the distribution to the data while keeping the ``loc``
-parameter null:
+The syntax ``distribution.fit(data, loc=1)`` can be used to fit the distribution to the data while keeping the ``loc``
+parameter at a fixed value, in this case 1:
 
->>> Normal.fit(data, loc=1)
+>>> Normal().fit(data, loc=1)
 Normal(loc=1.0, scale=1.9377929687500024)

docs/source/user_guide/penalty_fitting.rst

@@ -13,11 +13,11 @@ apply a Lasso penalty.
 >>> def lassolike_score(distribution, data):
 ...     return -np.sum(distribution.logpdf(data)) + 5 * np.abs(distribution.loc())
 ...
->>> cond_fit = Normal.fit(data, score=lassolike_score)
+>>> cond_fit = Normal().fit(data, score=lassolike_score)
 
 We then compare a fit using the standard negative log-likelihood function to the use of the Lasso-penalized likelihood.
 
->>> std_fit = Normal.fit(data)
+>>> std_fit = Normal().fit(data)
 >>> std_fit.loc.value
 -0.010891307380632494
 >>> cond_fit.loc.value

docs/source/user_guide/trend_fitting.rst

@@ -12,13 +12,13 @@ array([-0.99802364, -0.99503679, -0.98900434, -0.98277981, -0.979487  ,
 
 If we try to fit this without a trend, the resulting distribution will miss out on most of the information.
 
->>> Normal.fit(data)
+>>> Normal().fit(data)
 Normal(loc=-3.6462053656578005e-05, scale=0.5789668679237372)
 
 Fitting a ``Normal`` distribution with a trend in the ``loc`` parameter can be done using the following piece of code:
 
 >>> from pykelihood import kernels
->>> Normal.fit(data, loc=kernels.linear(np.arange(365)))
+>>> Normal().fit(data, loc=kernels.linear(np.arange(365)))
 Normal(loc=linear(a=-1.0000458359290572, b=0.005494714384381866), scale=0.0010055323717468906)
 
 The ``kernels`` module is flexible and can be adapted by users to support any kind of trend.

pykelihood/distributions/base.py

@@ -1,5 +1,6 @@
 from __future__ import annotations
 
+import warnings
 from abc import ABC, abstractmethod
 from collections.abc import Sequence
 from dataclasses import dataclass
@@ -12,7 +13,7 @@
 
 from pykelihood.generic_types import Obs
 from pykelihood.metrics import opposite_log_likelihood
-from pykelihood.parameters import ConstantParameter, Parametrized, ensure_parametrized
+from pykelihood.parameters import Parametrized, ensure_parametrized
 
 if TYPE_CHECKING:
     from typing import Self
@@ -49,8 +50,6 @@ class Distribution(Parametrized, ABC):
         Inverse of the cumulative distribution function.
     fit(data: Obs, x0: Sequence[float] = None, score: Callable[["Distribution", Obs], float] = opposite_log_likelihood, scipy_args: Optional[Dict] = None, **fixed_values) -> SomeDistribution
         Fit the distribution to the data.
-    fit_instance(data, score=opposite_log_likelihood, x0: Sequence[float] = None, scipy_args: Optional[Dict] = None, **fixed_values)
-        Fit the instance to the data.
     """
 
     def __hash__(self):
@@ -101,15 +100,14 @@ def inverse_cdf(self, q: Obs):
     def _apply_constraints(self, data):
         return data
 
-    @classmethod
     def fit(
-        cls: type[SomeDistribution],
+        self,
         data: Obs,
         x0: Sequence[float] | None = None,
         score: Callable[[Distribution, Obs], float] = opposite_log_likelihood,
         scipy_args: dict | None = None,
         **fixed_values,
-    ) -> Fit[SomeDistribution]:
+    ) -> Fit[Self]:
         """
         Fit the distribution to the data.
 
@@ -128,25 +126,14 @@ def fit(
 
         Returns
         -------
-        The result of the fit
+        The result of the fit. A new instance is created with the fitted parameters.
         """
-        init_parms = {}
-        for k in cls.params_names:
-            if k in fixed_values:
-                v = fixed_values.pop(k)
-                if isinstance(v, Parametrized):
-                    init_parms[k] = v
-                else:
-                    init_parms[k] = ConstantParameter(v)
-        # Add keyword arguments useful for object creation
-        for k, v in fixed_values.items():
-            if k not in init_parms:
-                init_parms[k] = v
-        init = cls(**init_parms)
+        init_parms = self._process_fit_params(**fixed_values)
+        init = type(self)(**init_parms)
         data = init._apply_constraints(data)
 
         if x0 is None:
-            x0 = [x.value for x in init.optimisation_params]
+            x0 = [x() for x in init.optimisation_params]
         else:
             if len(x0) != len(init.optimisation_params):
                 raise ValueError(
@@ -167,6 +154,14 @@ def to_minimize(x) -> float:
 
         return Fit(dist, data, score, x0=x0, optimize_result=optimization_result)
 
+    def fit_instance(self, *args, **kwargs):
+        warnings.warn(
+            "fit_instance is deprecated, use fit instead",
+            DeprecationWarning,
+            stacklevel=2,
+        )
+        return self.fit(*args, **kwargs)
+
     def _process_fit_params(self, **kwds):
         out_dict = self.param_dict.copy()
         to_remove = set()
@@ -195,38 +190,6 @@ def _process_fit_params(self, **kwds):
                 out_dict[name] = value
         return out_dict
 
-    def fit_instance(
-        self,
-        data: Obs,
-        score=opposite_log_likelihood,
-        x0: Sequence[float] | None = None,
-        scipy_args: dict | None = None,
-        **fixed_values,
-    ) -> Fit[Self]:
-        """
-        Fit the instance to the data.
-
-        Parameters
-        ----------
-        data : Obs
-            Data to fit the instance to.
-        score : Callable[["Distribution", Obs], float], optional
-            Scoring function, by default opposite_log_likelihood.
-        x0 : Sequence[float], optional
-            Initial guess for the parameters, by default None.
-        scipy_args : Optional[Dict], optional
-            Additional arguments for scipy.optimize.minimize, by default None.
-        fixed_values : dict
-            Fixed values for the parameters.
-
-        Returns
-        -------
-        Distribution
-            Fitted instance.
-        """
-        param_dict = self._process_fit_params(**fixed_values)
-        return self.fit(data, score=score, x0=x0, scipy_args=scipy_args, **param_dict)
-
 
 @dataclass
 class Fit(Generic[_T]):

pykelihood/distributions/custom.py

@@ -425,7 +425,7 @@ def _apply_constraints(self, x):
         """
         return x[self._valid_indices(x)]
 
-    def fit_instance(self, *args, **kwargs):
+    def fit(self, *args, **kwargs):
         """
         Fit the instance to the data.
 
@@ -442,7 +442,7 @@ def fit_instance(self, *args, **kwargs):
             The fitted instance.
         """
         kwargs.update(lower_bound=self.lower_bound, upper_bound=self.upper_bound)
-        return super().fit_instance(*args, **kwargs)
+        return super().fit(*args, **kwargs)
 
     def rvs(self, size: int, *args, **kwargs):
         """

pykelihood/parameters.py

@@ -166,14 +166,13 @@ def param_mapping(self, only_opt=False):
         return results
 
     @property
-    def optimisation_params(self) -> tuple[Parametrized]:
+    def optimisation_params(self) -> tuple[Parameter, ...]:
         """
         Get all parameters used in the optimization.
 
         Returns
         -------
-        Tuple[Parametrized]
-            The optimization parameters.
+        The optimization parameters.
         """
         unique = []
         for q in (p_ for p in self.params for p_ in p.optimisation_params):
@@ -182,13 +181,13 @@ def optimisation_params(self) -> tuple[Parametrized]:
         return unique
 
     @property
-    def optimisation_param_dict(self) -> dict[str, Parametrized]:
+    def optimisation_param_dict(self) -> dict[str, Parameter]:
         """
         Get a dictionary of optimization parameter names and their values.
 
         Returns
         -------
-        Dict[str, Parametrized]
+        Dict[str, Parameter]
             The optimization parameter dictionary.
         """
         p_dict = flatten_dict(self._optimisation_param_dict_helper())

pykelihood/profiler.py

@@ -86,7 +86,7 @@ def optimum(self):
         tuple
             A tuple containing the estimate and the score function value.
         """
-        estimate = self.distribution.fit_instance(
+        estimate = self.distribution.fit(
             self.data,
             score=self.score_function,
             x0=self.x0,
@@ -148,11 +148,7 @@ def test_profile_likelihood(self, range_for_param, param):
         profile_ll = []
         params = []
         for x in range_for_param:
-            pl = opt.fit_instance(
-                self.data,
-                score=self.score_function,
-                **{param: x},
-            )
+            pl = opt.fit(self.data, score=self.score_function, **{param: x})
             pl_value = -self.score_function(pl, self.data)
             pl_value = pl_value if isinstance(pl_value, float) else pl_value[0]
             if np.isfinite(pl_value):
@@ -190,9 +186,7 @@ def confidence_interval(self, param: str, precision=1e-5) -> tuple[float, float]
         value_threshold = func - chi2.ppf(self.inference_confidence, df=1) / 2
 
         def score(x: float):
-            new_opt = opt.fit_instance(
-                self.data, score=self.score_function, **{param: x}
-            )
+            new_opt = opt.fit(self.data, score=self.score_function, **{param: x})
             return -self.score_function(new_opt, self.data)
 
         def delta_to_threshold(x: float):

pyproject.toml

@@ -57,7 +57,3 @@ extend-select = [
     "I", # Import sorting
     "UP",  # PyUpgrade
 ]
-
-[tool.pytest.ini_options]
-durations = 10
-verbose = true

tests/test_distributions.py

@@ -29,21 +29,22 @@ class TestGEV:
     def test_fit(self, datasets):
         for ds in datasets:
             c, loc, scale = stats.genextreme.fit(ds)
-            fit = GEV.fit(ds)
+            fit = GEV().fit(ds)
             assert fit.loc() == approx(loc)
             assert fit.scale() == approx(scale)
             assert fit.shape() == approx(-c)
 
     def test_fixed_values(self):
         data = np.random.standard_normal(1000)
-        raw = Normal.fit(data)
+        raw = Normal().fit(data)
         assert raw.loc() == approx(0.0)
         assert raw.scale() == approx(1.0)
-        fixed = Normal.fit(data, loc=1.0)
+        fixed = Normal().fit(data, loc=1.0)
         assert fixed.loc() == 1.0
 
 
 def test_cache():
+    """There is no cache anymore, the test is kept as it can still be useful."""
     n = Normal(0, 1)
     np.testing.assert_array_almost_equal(
         n.pdf([-1, 0, 1]), [0.24197072, 0.39894228, 0.24197072]
@@ -87,28 +88,27 @@ def test_named_with_params_partial_assignment():
     assert m.scale() == 3
 
 
-def test_fit_instance(dataset):
-    std_fit = Normal.fit(dataset)
-    instance_fit = Normal(loc=kernels.constant()).fit_instance(dataset)
-    assert std_fit.loc() == approx(instance_fit.loc())
+def test_simple_fit(dataset):
+    std_fit = Normal().fit(dataset)
+    kernel_fit = Normal(loc=kernels.constant()).fit(dataset)
+    assert std_fit.loc() == approx(kernel_fit.loc())
 
 
-def test_fit_instance_fixed_params(dataset):
-    n = Normal().fit_instance(dataset, loc=5)
+def test_fit_fixed_param(dataset):
+    n = Normal().fit(dataset, loc=5)
     assert n.loc() == 5
 
 
-def test_fit_instance_fixed_params_multi_level(dataset, linear_kernel):
+def test_fit_fixed_param_depth_2(dataset, linear_kernel):
     n = Normal(loc=linear_kernel)
-    m = n.fit_instance(dataset, loc_a=5)
+    m = n.fit(dataset, loc_a=5)
     assert m.loc.a() == 5
 
 
-def test_fit_instance_fixed_params_extra_levels(dataset):
+def test_fit_fixed_param_depth_3(dataset):
     covariate = np.arange(len(dataset))
     n = Normal(loc=kernels.linear(covariate, a=kernels.linear(covariate)))
-    n.param_mapping()
-    m = n.fit_instance(dataset, loc_a_a=5)
+    m = n.fit(dataset, loc_a_a=5)
     assert m.loc.a.a() == 5
 
 
@@ -141,8 +141,8 @@ def test_truncated_distribution_fit():
     data = n.rvs(10000)
     trunc_data = data[data >= 0]
     truncated = TruncatedDistribution(Normal(), lower_bound=0)
-    fitted_all_data = truncated.fit_instance(data)
-    fitted_trunc = truncated.fit_instance(trunc_data)
+    fitted_all_data = truncated.fit(data)
+    fitted_trunc = truncated.fit(trunc_data)
     for p_trunc, p_all in zip(
         fitted_trunc.flattened_params, fitted_all_data.flattened_params
     ):
@@ -156,27 +156,27 @@ def test_distribution_fit_with_shared_params_in_trends():
     """
     x = np.array(np.random.uniform(size=200))
     y = np.array(np.random.normal(size=200))
-    alpha0_init = 0.0
-    alpha = Parameter(alpha0_init)
-    n = Normal.fit(y, loc=linear(x=x, b=alpha), scale=linear(x=x, b=alpha))
+    alpha = Parameter(0.0)
+    n = Normal().fit(y, loc=linear(x=x, b=alpha), scale=linear(x=x, b=alpha))
     alpha1 = n.loc.b
     alpha2 = n.scale.b
     assert alpha1 == alpha2
 
 
-def test_fit_instance_fixing_shared_params_in_trends():
+def test_fit_fixing_shared_params_in_trends():
     """
-    when 2 trends in the distribution parameters share a common parameter, e.g. alpha in the below example, making one of the corresponding trend parameter constant should automatically result in the other trend parameter is constant.
+    when 2 trends in the distribution parameters share a common parameter,
+    e.g. alpha in the below example, making one of the corresponding trend parameter
+    constant should automatically result in the other trend parameter being constant.
     """
     x = np.array(np.random.uniform(size=200))
     y = np.array(np.random.normal(size=200))
-    alpha0_init = 0.0
-    alpha = Parameter(alpha0_init)
-    n = Normal.fit(y, loc=linear(x=x, b=alpha), scale=linear(x=x, b=alpha))
+    alpha = Parameter(0.0)
+    n = Normal().fit(y, loc=linear(x=x, b=alpha), scale=linear(x=x, b=alpha))
     fixed_alpha = ConstantParameter(
         n.loc.b.value
-    )  # should be equal to fit.scale.a as per problem1
-    fit_with_fixed_alpha = n.fit_instance(data=y, loc_b=fixed_alpha)
+    )  # should be equal to fit.scale.b as per previous test above
+    fit_with_fixed_alpha = n.fit(data=y, loc_b=fixed_alpha)
     assert isinstance(fit_with_fixed_alpha.scale.b, ConstantParameter)
     assert fit_with_fixed_alpha.scale.b.value == fixed_alpha.value