Addition of non-uniform priors.

docs/source/arrhenius_t.ipynb

@@ -31,6 +31,7 @@
     "import matplotlib.pyplot as plt\n",
     "from kinisi.analyze import DiffusionAnalyzer\n",
     "from kinisi.arrhenius import Arrhenius\n",
+    "from scipy.stats import uniform\n",
     "import warnings\n",
     "np.random.seed(42)\n",
     "warnings.filterwarnings(\"ignore\", category=UserWarning)\n",
@@ -145,26 +146,31 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "s = Arrhenius(td, bounds=[[0.1 * sc.Unit('eV'), 0.2 * sc.Unit('eV')], \n",
-    "                          [1e-5 * td.data.unit, 1e-4 * td.data.unit]])"
+    "priors = [uniform(0.1, 0.1), uniform(1e-5, 9e-5)]\n",
+    "\n",
+    "s = Arrhenius(td, priors=priors)"
    ]
   },
   {
    "cell_type": "markdown",
    "id": "507098fb-3afe-40b7-977d-1db0d8eb20a4",
    "metadata": {},
    "source": [
-    "When you create a `TemperatureDependent` object, of which `Arrhenius` is a subclass, you can optionally provide bounds as a list of lists. \n",
+    "When you create a `TemperatureDependent` object, of which `Arrhenius` is a subclass, you can optionally provide priors as a list of `scipy.stats` objects. \n",
     "For example, if we wanted the activation energy to be restricted between 0.1 and 0.2 eV and the preexponential factor between 10<sup>-5</sup> and 10<sup>-4</sup> cm<sup>2</sup>/s, then we would use the following: \n",
     "\n",
     "```python\n",
-    "Arrhenius(td, bounds=[[0.1 * sc.Unit('eV'), 0.2 * sc.Unit('eV')], \n",
-    "                      [1e-5 * td.data.unit, 1e-4 * td.data.unit]])\n",
-    "```\n",
+    "priors = [uniform(0.1, 0.1), uniform(1e-5, 9e-5)]\n",
     "\n",
-    "If you do not set bounds, they are automatically set to 0.5 and 1.5 of the best fit values.\n",
+    "s = Arrhenius(td, priors=priors)\n",
+    "```\n",
     "\n",
-    "The `Arrhenius` object automatically estimates the maximum likelihood values for the activation energy and preexponential factor and stores them is an appropriate `sc.DataGroup`. "
+    "So there are two uniform distributions, one for the activation energy and another for the preexponential factor.\n",
+     "The first value in `uniform` is `loc` (the minimum), and the second is `scale` (the width, i.e., max - min).  \n",
+     "We recommend looking at the `scipy.stats` documentation to clarify how to use these objects. \n",
+     "If you do not set `priors`, they are automatically set to be uniform with a min and max of 0.5 and 1.5 times the maximum-likelihood values.\n",
+     "\n",
+     "The `Arrhenius` object automatically estimates the maximum-likelihood values for the activation energy and preexponential factor and stores them in an appropriate `sc.DataGroup`. "
    ]
   },
   {
@@ -308,7 +314,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.13"
+   "version": "3.13.11"
   }
  },
  "nbformat": 4,

kinisi/__init__.py

@@ -1,4 +1,4 @@
-__version__ = '2.0.5'
+__version__ = '2.1.0'
 
 from .due import Doi, due
 

kinisi/arrhenius.py

@@ -33,8 +33,9 @@ class TemperatureDependent(FittingBase):
     :param function: A callable function that describes the relationship between temperature and diffusion.
     :param parameter_names: A tuple of parameter names for the function.
     :param parameter_units: A tuple of sc.Unit objects corresponding to the parameter names.
-    :param bounds: Optional bounds for the parameters of the function. Defaults to None, in which case these
-        are defined as +/- 50 percent of the best fit values.
+    :param priors: Optional prior probability distributions for the parameters of the function.
+        Defaults to None, in which case a uniform distribution is defined with limits of
+        +/- 50 percent of the best fit values.
     """
 
     def __init__(
@@ -43,7 +44,7 @@ def __init__(
         function: Callable,
         parameter_names: tuple[str],
         parameter_units: tuple[sc.Unit],
-        bounds: None | list = None,
+        priors: None | list = None,
     ) -> 'TemperatureDependent':
         self.diffusion = diffusion
         self.temperature = diffusion.coords['temperature']
@@ -53,7 +54,7 @@ def __init__(
             function=function,
             parameter_names=parameter_names,
             parameter_units=parameter_units,
-            bounds=bounds,
+            priors=priors,
             coordinate_name='temperature',
         )
 
@@ -94,19 +95,20 @@ class Arrhenius(TemperatureDependent):
     Evaluate the data with a standard Arrhenius relationship.
 
     :param diffusion: Diffusion coefficient sc.DataFrame with a temperature coordinate and variances.
-    :param bounds: Optional bounds for the parameters of the function. Defaults to None, in which case these
-        are defined as +/- 50 percent of the best fit values.
+    :param priors: Optional prior probability distributions for the parameters of the function.
+        Defaults to None, in which case a uniform distribution is defined with limits of
+        +/- 50 percent of the max likelihood values.
     """
 
     def __init__(
         self,
         diffusion,
-        bounds: tuple[tuple[VariableLike, VariableLike], tuple[VariableLike, VariableLike]] | None = None,
+        priors: list | None = None,
     ) -> 'Arrhenius':
         parameter_names = ('activation_energy', 'preexponential_factor')
         parameter_units = (sc.Unit('eV'), sc.Unit('cm^2/s'))
 
-        super().__init__(diffusion, arrhenius, parameter_names, parameter_units, bounds=bounds)
+        super().__init__(diffusion, arrhenius, parameter_names, parameter_units, priors=priors)
 
     @property
     def activation_energy(self) -> VariableLike | Samples:
@@ -144,22 +146,20 @@ class VogelFulcherTammann(TemperatureDependent):
     Evaluate the data with a Vogel-Fulcher-Tammann relationship.
 
     :param diffusion: Diffusion coefficient sc.DataFrame with a temperature coordinate and variances.
-    :param bounds: Optional bounds for the parameters of the function. Defaults to None, in which case these
-        are defined as +/- 50 percent of the best fit values.
+    :param priors: Optional prior probability distributions for the parameters of the function.
+        Defaults to None, in which case a uniform distribution is defined with limits of
+        +/- 50 percent of the max likelihood values.
     """
 
     def __init__(
         self,
         diffusion,
-        bounds: tuple[
-            tuple[VariableLike, VariableLike], tuple[VariableLike, VariableLike], tuple[VariableLike, VariableLike]
-        ]
-        | None = None,
+        priors: list | None = None,
     ) -> 'VogelFulcherTammann':
         parameter_names = ('activation_energy', 'preexponential_factor', 'T0')
         parameter_units = (sc.Unit('eV'), sc.Unit('cm^2/s'), sc.Unit('K'))
 
-        super().__init__(diffusion, vtf_equation, parameter_names, parameter_units, bounds=bounds)
+        super().__init__(diffusion, vtf_equation, parameter_names, parameter_units, priors=priors)
 
     @property
     def activation_energy(self) -> VariableLike | Samples:

kinisi/fitting.py

@@ -31,8 +31,9 @@ class FittingBase:
     :param function: A callable function that describes the relationship
     :param parameter_names: A tuple of parameter names for the function
     :param parameter_units: A tuple of sc.Unit objects corresponding to the parameter names
-    :param bounds: Optional bounds for the parameters of the function. Defaults to None,
-        in which case these are defined as +/- 50 percent of the best fit values
+    :param priors: Prior probability distributions for the parameters of the function.
+        Defaults to None, in which case a uniform distribution is defined with limits of
+        +/- 50 percent of the max likelihood fit values.
     :param coordinate_name: Name of the coordinate to use as independent variable
     """
 
@@ -42,28 +43,36 @@ def __init__(
         function: Callable,
         parameter_names: tuple[str, ...],
         parameter_units: tuple[sc.Unit, ...],
-        bounds: None | list = None,
+        priors: None | list = None,
         coordinate_name: str | None = None,
     ) -> 'FittingBase':
         self.data = data
         self.data_group = sc.DataGroup({'data': data})
         self.function = function
-        self.bounds = bounds
+        self.priors = priors
         self.parameter_names = parameter_names
         self.parameter_units = parameter_units
         self.coordinate_name = coordinate_name
 
-        if bounds is not None and len(bounds) != len(self.parameter_names):
-            raise ValueError(
-                f'Bounds must be a tuple of length {len(self.parameter_names)}, got {len(bounds)} instead.'
-            )
-
-        # Perform initial fit
-        self.max_likelihood()
+        if priors is not None:
+            if len(priors) != len(self.parameter_names):
+                raise ValueError(
+                    f'Priors must be a list of length {len(self.parameter_names)}, got {len(priors)} instead.'
+                )
+            if not all(hasattr(p, 'logpdf') and hasattr(p, 'ppf') for p in priors):
+                raise ValueError(
+                    "Priors must provide 'logpdf' and 'ppf' methods (e.g., frozen scipy.stats distributions)."
+                )
+
+        if self.priors is None:
+            # Perform initial fit
+            self.max_likelihood()
+        else:
+            self.max_aposteriori()
 
-        # Set default bounds if not provided
-        if self.bounds is None:
-            self.bounds = tuple(
+        # Set default priors if not provided
+        if self.priors is None:
+            bounds = tuple(
                 [
                     (
                         self.data_group[p].value * 0.5 * self.data_group[p].unit,
@@ -73,8 +82,8 @@ def __init__(
                 ]
             )
 
-        # Set up priors for nested sampling
-        self.priors = [uniform(b[0].value, b[1].value - b[0].value) for b in self.bounds]
+            # Set up priors for nested sampling
+            self.priors = [uniform(b[0].value, b[1].value - b[0].value) for b in bounds]
 
     def __repr__(self):
         """String representation."""
@@ -151,6 +160,15 @@ def log_posterior(self, parameters: tuple[float]) -> float:
         """
         return self.log_likelihood(parameters) + self.log_prior(parameters)
 
+    def nlp(self, parameters: tuple[float]) -> float:
+        """
+        Calculate the negative log posterior of the model given the data.
+
+        :param parameters: The parameters of the model.
+        :return: The negative log posterior of the model.
+        """
+        return -self.log_posterior(parameters)
+
     def prior_transform(self, parameters: tuple[float]) -> tuple[float]:
         """
         Transform the parameters from the prior space to the parameter space.
@@ -163,14 +181,23 @@ def prior_transform(self, parameters: tuple[float]) -> tuple[float]:
             x[i] = self.priors[i].ppf(parameters[i])
         return x
 
-    def max_likelihood(self) -> tuple[float]:
-        """Find the best fit parameters for the model."""
-        if self.bounds is not None:
-            x0 = [((b[1] + b[0]) / 2) for b in self.bounds]
+    def max_likelihood(self) -> None:
+        """Find the max likelihood fit parameters for the model."""
+        if self.priors is not None:
+            x0 = [p.mean() for p in self.priors]
+        else:
+            x0 = [1 for u in self.parameter_units]
+        result = minimize(self.nll, x0).x
+        for i, name in enumerate(self.parameter_names):
+            self.data_group[name] = result[i] * self.parameter_units[i]
+
+    def max_aposteriori(self) -> None:
+        """Find the max aposteriori fit parameters for the model."""
+        if self.priors is not None:
+            x0 = [p.mean() for p in self.priors]
         else:
-            x0 = [1 * u for u in self.parameter_units]
-        bounds = [(b[0].value, b[1].value) for b in self.bounds] if self.bounds is not None else None
-        result = minimize(self.nll, [x.value for x in x0], bounds=bounds).x
+            x0 = [1 for u in self.parameter_units]
+        result = minimize(self.nlp, x0).x
         for i, name in enumerate(self.parameter_names):
             self.data_group[name] = result[i] * self.parameter_units[i]