diff --git a/docs/source/arrhenius_t.ipynb b/docs/source/arrhenius_t.ipynb
index 74f7e671..1564ddb9 100644
--- a/docs/source/arrhenius_t.ipynb
+++ b/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,8 +146,9 @@
"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)"
]
},
{
@@ -154,17 +156,21 @@
"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-5 and 10-4 cm2/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,
diff --git a/kinisi/__init__.py b/kinisi/__init__.py
index 5683f0eb..99358322 100644
--- a/kinisi/__init__.py
+++ b/kinisi/__init__.py
@@ -1,4 +1,4 @@
-__version__ = '2.0.5'
+__version__ = '2.1.0'
from .due import Doi, due
diff --git a/kinisi/arrhenius.py b/kinisi/arrhenius.py
index 6291e2a4..7c21dda7 100644
--- a/kinisi/arrhenius.py
+++ b/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:
diff --git a/kinisi/fitting.py b/kinisi/fitting.py
index 5ff15737..82d5e4b1 100644
--- a/kinisi/fitting.py
+++ b/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]
diff --git a/kinisi/tests/test_arrhenius.py b/kinisi/tests/test_arrhenius.py
index b5518141..32aef08d 100644
--- a/kinisi/tests/test_arrhenius.py
+++ b/kinisi/tests/test_arrhenius.py
@@ -10,10 +10,8 @@
import unittest
-import numpy as np
import scipp as sc
from numpy.testing import assert_almost_equal, assert_equal
-from scipp import testing
from scipy.stats import uniform
from kinisi import arrhenius
@@ -39,129 +37,9 @@ def straight_line(x, m, c):
class TestTemperatureDependent(unittest.TestCase):
"""
- Unit tests for TemperatureDependent class
+ Unit tests for the TemperatureDependent class
"""
- def test_init(self):
- """
- Test the initialisation of TemperatureDependent class
- """
- td = arrhenius.TemperatureDependent(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')))
- u = uniform(loc=0, scale=1)
- assert_equal(td.function, straight_line)
- assert td.parameter_names == ('m', 'c')
- assert td.parameter_units == (sc.Unit('m/s'), sc.Unit('m'))
- assert isinstance(td.data_group['m'], sc.Variable)
- assert isinstance(td.data_group['c'], sc.Variable)
- assert isinstance(td.bounds[0][0], sc.Variable)
- assert isinstance(td.bounds[0][1], sc.Variable)
- assert isinstance(td.bounds[1][0], sc.Variable)
- assert isinstance(td.bounds[1][1], sc.Variable)
- assert isinstance(td.priors[0], type(u))
- assert isinstance(td.priors[1], type(u))
-
- def test_init_bounds(self):
- """
- Test the initialisation of TemperatureDependent class with bounds
- """
- bounds = ((0 * sc.Unit('m/s'), 1 * sc.Unit('m/s')), (0 * sc.Unit('m'), 1e20 * sc.Unit('m')))
- td = arrhenius.TemperatureDependent(
- data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')), bounds=bounds
- )
- assert_equal(td.function, straight_line)
- assert td.parameter_names == ('m', 'c')
- assert td.parameter_units == (sc.Unit('m/s'), sc.Unit('m'))
- testing.assert_allclose(td.bounds[0][0], bounds[0][0])
- testing.assert_allclose(td.bounds[0][1], bounds[0][1])
- testing.assert_allclose(td.bounds[1][0], bounds[1][0])
- testing.assert_allclose(td.bounds[1][1], bounds[1][1])
-
- def test_init_wrong_bounds(self):
- """
- Test the initialisation of TemperatureDependent class with wrong bounds
- """
- with self.assertRaises(ValueError):
- arrhenius.TemperatureDependent(
- data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')), bounds=((0, 1), (0, 1e20), (0, 1))
- )
-
- def test_repr(self):
- """
- Test the string representation of TemperatureDependent class
- """
- td = arrhenius.TemperatureDependent(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')))
- assert str(td.__repr__()) == str(td.data_group.__repr__())
-
- def test_str(self):
- """
- Test the string representation of TemperatureDependent class
- """
- td = arrhenius.TemperatureDependent(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')))
- assert str(td) == str(td.data_group)
-
- def test_repr_html(self):
- """
- Test the HTML representation of TemperatureDependent class
- """
- td = arrhenius.TemperatureDependent(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')))
- assert isinstance(td._repr_html_(), str)
-
- def test_log_likelihood(self):
- """
- Test the log-likelihood function of TemperatureDependent class
- """
- td = arrhenius.TemperatureDependent(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')))
- assert isinstance(td.log_likelihood([1, 0]), float)
- assert_almost_equal(td.log_likelihood([1, 0]), -13.275855715784758)
-
- def test_nll(self):
- """
- Test the negative log-likelihood function of TemperatureDependent class
- """
- td = arrhenius.TemperatureDependent(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')))
- assert isinstance(td.nll([1, 0]), float)
- assert_almost_equal(td.nll([1, 0]), 13.275855715784758)
-
- def test_log_prior(self):
- """
- Test the log-prior function of TemperatureDependent class
- """
- td = arrhenius.TemperatureDependent(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')))
- assert isinstance(td.log_prior([1, 0]), float)
- assert_almost_equal(td.log_prior([1, 0]), -np.inf)
-
- def test_mcmc(self):
- """
- Test the MCMC sampling function of TemperatureDependent class
- """
- td = arrhenius.TemperatureDependent(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')))
- td.mcmc(n_samples=10, n_burn=5, n_walkers=32)
- assert isinstance(td.data_group['m'], Samples)
- assert isinstance(td.data_group['c'], Samples)
- assert td.data_group['m'].shape == (32,)
- assert td.data_group['c'].shape == (32,)
-
- def test_nested_sampling(self):
- """
- Test the nested sampling function of TemperatureDependent class
- """
- td = arrhenius.TemperatureDependent(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')))
- td.nested_sampling()
- assert isinstance(td.data_group['m'], Samples)
- assert isinstance(td.data_group['c'], Samples)
- assert isinstance(td.logz, sc.Variable)
-
- def test_flatchain(self):
- """
- Test the flatchain property of TemperatureDependent class
- """
- td = arrhenius.TemperatureDependent(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')))
- td.mcmc(n_samples=10, n_burn=5, n_walkers=32)
- td.mcmc(n_samples=10, n_burn=5, n_walkers=32)
- assert isinstance(td.flatchain, sc.DataGroup)
- assert len(td.flatchain) == 2
- assert td.flatchain.shape == (32,)
-
def test_extrapolate(self):
"""
Test the extrapolate function of TemperatureDependent class
@@ -199,19 +77,17 @@ def test_init(self):
assert isinstance(arr.activation_energy, sc.Variable)
assert isinstance(arr.preexponential_factor, sc.Variable)
- def test_init_bounds(self):
+ def test_init_priors(self):
"""
- Test the initialisation of Arrhenius class with bounds
+ Test the initialisation of Arrhenius class with priors
"""
- bounds = ((0 * sc.Unit('eV'), 1 * sc.Unit('eV')), (0 * sc.Unit('cm^2/s'), 1e20 * sc.Unit('cm^2/s')))
- arr = arrhenius.Arrhenius(data, bounds=bounds)
+ priors = [uniform(0, 1), uniform(0, 1e20)]
+ arr = arrhenius.Arrhenius(data, priors=priors)
assert_equal(arr.function, arrhenius.arrhenius)
assert arr.parameter_names == ('activation_energy', 'preexponential_factor')
assert arr.parameter_units == (sc.Unit('eV'), sc.Unit('cm^2/s'))
- testing.assert_allclose(arr.bounds[0][0], bounds[0][0])
- testing.assert_allclose(arr.bounds[0][1], bounds[0][1])
- testing.assert_allclose(arr.bounds[1][0], bounds[1][0])
- testing.assert_allclose(arr.bounds[1][1], bounds[1][1])
+ assert arr.priors[0] == priors[0]
+ assert arr.priors[1] == priors[1]
def test_arrhenius(self):
"""
@@ -237,25 +113,22 @@ def test_init(self):
assert isinstance(vtf.preexponential_factor, sc.Variable)
assert isinstance(vtf.T0, sc.Variable)
- def test_init_bounds(self):
+ def test_init_priors(self):
"""
- Test the initialisation of VogelFulcherTammann class with bounds
+ Test the initialisation of VogelFulcherTammann class with priors
"""
- bounds = (
- (0 * sc.Unit('eV'), 1 * sc.Unit('eV')),
- (0 * sc.Unit('cm^2/s'), 1e20 * sc.Unit('cm^2/s')),
- (0 * sc.Unit('K'), 1000 * sc.Unit('K')),
+ priors = (
+ uniform(0, 1),
+ uniform(0, 1e20),
+ uniform(0, 1000),
)
- vtf = arrhenius.VogelFulcherTammann(data, bounds=bounds)
+ vtf = arrhenius.VogelFulcherTammann(data, priors=priors)
assert_equal(vtf.function, arrhenius.vtf_equation)
assert vtf.parameter_names == ('activation_energy', 'preexponential_factor', 'T0')
assert vtf.parameter_units == (sc.Unit('eV'), sc.Unit('cm^2/s'), sc.Unit('K'))
- testing.assert_allclose(vtf.bounds[0][0], bounds[0][0])
- testing.assert_allclose(vtf.bounds[0][1], bounds[0][1])
- testing.assert_allclose(vtf.bounds[1][0], bounds[1][0])
- testing.assert_allclose(vtf.bounds[1][1], bounds[1][1])
- testing.assert_allclose(vtf.bounds[2][0], bounds[2][0])
- testing.assert_allclose(vtf.bounds[2][1], bounds[2][1])
+ assert vtf.priors[0] == priors[0]
+ assert vtf.priors[1] == priors[1]
+ assert vtf.priors[2] == priors[2]
def test_vtf_equation(self):
"""
diff --git a/kinisi/tests/test_fitting.py b/kinisi/tests/test_fitting.py
new file mode 100644
index 00000000..322c91f8
--- /dev/null
+++ b/kinisi/tests/test_fitting.py
@@ -0,0 +1,177 @@
+"""
+Tests for fitting module
+"""
+
+# Copyright (c) kinisi developers.
+# Distributed under the terms of the MIT License.
+# author: Andrew R. McCluskey (arm61)
+
+# pylint: disable=R0201
+
+import unittest
+
+import numpy as np
+import scipp as sc
+from numpy.testing import assert_almost_equal, assert_equal
+from scipy.stats import norm, uniform
+from scipy.stats._distn_infrastructure import rv_frozen
+
+from kinisi import fitting
+from kinisi.samples import Samples
+
+temp = sc.linspace(start=5, stop=50, num=10, dim='temperature', unit=sc.units.K)
+D = sc.linspace(start=5, stop=50, num=10, unit='cm^2 / s', dim='temperature')
+D.variances = D.values * 0.1 # 10% uncertainty
+data = sc.DataArray(data=D, coords={'temperature': temp})
+
+
+def straight_line(x, m, c):
+ """
+ A simple linear function for testing purposes.
+
+ :param x: The independent variable.
+ :param m: The slope of the line.
+ :param c: The y-intercept of the line.
+ :return: The value of the linear function at x.
+ """
+ return m * x + c
+
+
+class TestFittingBase(unittest.TestCase):
+ """
+ Unit tests for FittingBase class
+ """
+
+ def test_init(self):
+ """
+ Test the initialisation of FittingBase class
+ """
+ td = fitting.FittingBase(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')))
+ assert_equal(td.function, straight_line)
+ assert td.parameter_names == ('m', 'c')
+ assert td.parameter_units == (sc.Unit('m/s'), sc.Unit('m'))
+ assert isinstance(td.data_group['m'], sc.Variable)
+ assert isinstance(td.data_group['c'], sc.Variable)
+ assert isinstance(td.priors[0], rv_frozen)
+ assert isinstance(td.priors[1], rv_frozen)
+
+ def test_init_priors(self):
+ """
+ Test the initialisation of FittingBase class with priors
+ """
+ priors = [uniform(0, 1), uniform(0, 1e20)]
+ td = fitting.FittingBase(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')), priors=priors)
+ assert_equal(td.function, straight_line)
+ assert td.parameter_names == ('m', 'c')
+ assert td.parameter_units == (sc.Unit('m/s'), sc.Unit('m'))
+ assert td.priors[0] == priors[0]
+ assert td.priors[1] == priors[1]
+
+ def test_init_norm_priors(self):
+ """
+ Test the initialisation of the FittingBase Class
+ """
+ priors = [norm(10, 1), norm(0, 1)]
+ td = fitting.FittingBase(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')), priors=priors)
+ assert_equal(td.function, straight_line)
+ assert td.parameter_names == ('m', 'c')
+ assert td.parameter_units == (sc.Unit('m/s'), sc.Unit('m'))
+ assert td.priors[0] == priors[0]
+ assert td.priors[1] == priors[1]
+
+ def test_init_wrong_priors(self):
+ """
+ Test the initialisation of FittingBase class with wrong number of priors
+ """
+ with self.assertRaises(ValueError):
+ priors = [uniform(0, 1), uniform(0, 1e20), uniform(0, 1)]
+ fitting.FittingBase(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')), priors=priors)
+
+ def test_repr(self):
+ """
+ Test the string representation of FittingBase class
+ """
+ td = fitting.FittingBase(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')))
+ assert str(td.__repr__()) == str(td.data_group.__repr__())
+
+ def test_str(self):
+ """
+ Test the string representation of FittingBase class
+ """
+ td = fitting.FittingBase(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')))
+ assert str(td) == str(td.data_group)
+
+ def test_repr_html(self):
+ """
+ Test the HTML representation of FittingBase class
+ """
+ td = fitting.FittingBase(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')))
+ assert isinstance(td._repr_html_(), str)
+
+ def test_log_likelihood(self):
+ """
+ Test the log-likelihood function of FittingBase class
+ """
+ td = fitting.FittingBase(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')))
+ assert isinstance(td.log_likelihood([1, 0]), float)
+ assert_almost_equal(td.log_likelihood([1, 0]), -13.275855715784758)
+
+ def test_nll(self):
+ """
+ Test the negative log-likelihood function of FittingBase class
+ """
+ td = fitting.FittingBase(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')))
+ assert isinstance(td.nll([1, 0]), float)
+ assert_almost_equal(td.nll([1, 0]), 13.275855715784758)
+
+ def test_log_prior(self):
+ """
+ Test the log-prior function of FittingBase class
+ """
+ td = fitting.FittingBase(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')))
+ assert isinstance(td.log_prior([1, 0]), float)
+ assert_almost_equal(td.log_prior([1, 0]), -np.inf)
+
+ def test_mcmc(self):
+ """
+ Test the MCMC sampling function of FittingBase class
+ """
+ td = fitting.FittingBase(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')))
+ td.mcmc(n_samples=10, n_burn=5, n_walkers=32)
+ assert isinstance(td.data_group['m'], Samples)
+ assert isinstance(td.data_group['c'], Samples)
+ assert td.data_group['m'].shape == (32,)
+ assert td.data_group['c'].shape == (32,)
+
+ def test_nested_sampling(self):
+ """
+ Test the nested sampling function of FittingBase class
+ """
+ td = fitting.FittingBase(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')))
+ td.nested_sampling()
+ assert isinstance(td.data_group['m'], Samples)
+ assert isinstance(td.data_group['c'], Samples)
+ assert isinstance(td.logz, sc.Variable)
+
+ def test_nested_sampling_norm_priors(self):
+ """
+ Test the nested sampling function of FittingBase class where priors are norm
+ """
+ td = fitting.FittingBase(
+ data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')), priors=[norm(1, 1), norm(0, 1)]
+ )
+ td.nested_sampling()
+ assert isinstance(td.data_group['m'], Samples)
+ assert isinstance(td.data_group['c'], Samples)
+ assert isinstance(td.logz, sc.Variable)
+
+ def test_flatchain(self):
+ """
+ Test the flatchain property of FittingBase class
+ """
+ td = fitting.FittingBase(data, straight_line, ('m', 'c'), (sc.Unit('m/s'), sc.Unit('m')))
+ td.mcmc(n_samples=10, n_burn=5, n_walkers=32)
+ td.mcmc(n_samples=10, n_burn=5, n_walkers=32)
+ assert isinstance(td.flatchain, sc.DataGroup)
+ assert len(td.flatchain) == 2
+ assert td.flatchain.shape == (32,)
diff --git a/kinisi/tests/test_yeh_hummer.py b/kinisi/tests/test_yeh_hummer.py
index 92d48cf5..82159ae4 100644
--- a/kinisi/tests/test_yeh_hummer.py
+++ b/kinisi/tests/test_yeh_hummer.py
@@ -5,6 +5,7 @@
import numpy as np
import pytest
import scipp as sc
+from scipy.stats import uniform
from kinisi.yeh_hummer import YehHummer
@@ -91,8 +92,8 @@ def test_yeh_hummer_mcmc(self):
assert dist.shape[0] == len(box_lengths) # Number of data points
assert dist.shape[1] > 0 # Number of samples
- def test_yeh_hummer_bounds(self):
- """Test YehHummer with custom bounds."""
+ def test_yeh_hummer_priors(self):
+ """Test YehHummer with custom priors."""
box_lengths = np.array([20.0, 30.0, 40.0])
D_values = np.array([5.0e-5, 5.2e-5, 5.4e-5])
D_errors = np.array([0.1e-5, 0.1e-5, 0.1e-5])
@@ -102,17 +103,17 @@ def test_yeh_hummer_bounds(self):
coords={'box_length': sc.Variable(dims=['system'], values=box_lengths, unit='angstrom')},
)
- # Custom bounds
- bounds = (
- (4e-5 * sc.Unit('cm^2/s'), 7e-5 * sc.Unit('cm^2/s')), # D_0 bounds
- (1e-4 * sc.Unit('Pa*s'), 1e-2 * sc.Unit('Pa*s')), # viscosity bounds
+ # Custom priors
+ priors = (
+ uniform(4e-5, 7e-5), # D_0 prior
+ uniform(1e-4, 1e-2), # slope prior
)
- yh = YehHummer(td, temperature=sc.scalar(298, unit='K'), bounds=bounds)
+ yh = YehHummer(td, temperature=sc.scalar(298, unit='K'), priors=priors)
- # Check that fitted values are within bounds
- assert bounds[0][0].value <= yh.D_infinite.value <= bounds[0][1].value
- assert bounds[1][0].value <= yh.shear_viscosity.value <= bounds[1][1].value
+ # Check that fitted values are within priors
+ assert priors[0].a <= yh.D_infinite.value <= (priors[0].b + priors[0].a)
+ assert priors[1].a <= yh.shear_viscosity.value <= (priors[1].b + priors[1].a)
def test_yeh_hummer_properties(self):
"""Test YehHummer property accessors."""
@@ -137,8 +138,8 @@ def test_yeh_hummer_properties(self):
assert len(str(yh)) > 0
assert len(repr(yh)) > 0
- def test_yeh_hummer_invalid_bounds(self):
- """Test YehHummer with invalid bounds."""
+ def test_yeh_hummer_invalid_priors(self):
+ """Test YehHummer with invalid priors."""
box_lengths = np.array([20.0, 30.0, 40.0])
D_values = np.array([5.0e-5, 5.2e-5, 5.4e-5])
D_errors = np.array([0.1e-5, 0.1e-5, 0.1e-5])
@@ -149,7 +150,7 @@ def test_yeh_hummer_invalid_bounds(self):
)
# Wrong number of bounds
- bounds = ((4e-5 * sc.Unit('cm^2/s'), 7e-5 * sc.Unit('cm^2/s')),) # Only one bound
+ priors = [uniform(4e-5, 7e-5)] # Only one bound
- with pytest.raises(ValueError, match='Bounds must be a tuple of length 2'):
- YehHummer(td, temperature=sc.scalar(298, unit='K'), bounds=bounds)
+ with pytest.raises(ValueError, match='Priors must be a list of length 2.'):
+ YehHummer(td, temperature=sc.scalar(298, unit='K'), priors=priors)
diff --git a/kinisi/yeh_hummer.py b/kinisi/yeh_hummer.py
index ba2136e0..cb9b0113 100644
--- a/kinisi/yeh_hummer.py
+++ b/kinisi/yeh_hummer.py
@@ -13,6 +13,7 @@
import scipp as sc
import scipp.constants as const
from scipy.optimize import curve_fit, minimize
+from scipy.stats import uniform
from kinisi import __version__
from kinisi.fitting import FittingBase
@@ -40,10 +41,11 @@ class YehHummer(FittingBase):
:param diffusion: sc.DataArray with diffusion coefficients and box_length coordinate
:param temperature: Temperature (will be extracted from coords if not provided)
- :param bounds: Optional bounds for [D_0, viscosity] parameters (viscosity in Pa*s)
+ :param priors: Optional priors for [D_0, viscosity] parameters using scipy.stats objects,
+ where viscosity has units Pa*s.
"""
- def __init__(self, diffusion, temperature: sc.Variable, bounds=None):
+ def __init__(self, diffusion, temperature: sc.Variable, priors=None):
self.diffusion = diffusion
# Extract box lengths from coordinates
@@ -60,22 +62,22 @@ def __init__(self, diffusion, temperature: sc.Variable, bounds=None):
parameter_names = ('D_0', 'slope')
parameter_units = (diffusion.unit, self._slope_unit)
- # Compute bounds: use provided or defaults
- if bounds is None:
+ # Compute priors: use provided or defaults
+ if priors is None:
D_max = np.max(diffusion.values)
- D_bounds = (D_max * 0.8 * diffusion.unit, D_max * 2.0 * diffusion.unit)
+ D_prior = uniform(D_max * 0.8, D_max * 1.2)
visc_lower, visc_upper = 1e-5 * sc.Unit('Pa*s'), 1e-1 * sc.Unit('Pa*s')
+
+ # Higher viscosity = lower slope, so bounds are inverted
+ slope_bounds = (
+ self.viscosity_to_slope(visc_upper) * self._slope_unit,
+ self.viscosity_to_slope(visc_lower) * self._slope_unit,
+ )
+ slope_prior = uniform(slope_bounds[0].value, (slope_bounds[1] - slope_bounds[0]).value)
+ priors = [D_prior, slope_prior]
else:
- if len(bounds) != 2:
- raise ValueError('Bounds must be a tuple of length 2: (D_0_bounds, viscosity_bounds)')
- D_bounds = (bounds[0][0].to(unit=parameter_units[0]), bounds[0][1].to(unit=parameter_units[0]))
- visc_lower, visc_upper = bounds[1]
-
- # Higher viscosity = lower slope, so bounds are inverted
- slope_bounds = (
- self.viscosity_to_slope(visc_upper) * self._slope_unit,
- self.viscosity_to_slope(visc_lower) * self._slope_unit,
- )
+ if len(priors) != 2:
+ raise ValueError('Priors must be a list of length 2.')
# Initialize base class with linear function
super().__init__(
@@ -83,7 +85,7 @@ def __init__(self, diffusion, temperature: sc.Variable, bounds=None):
function=self._yeh_hummer_function,
parameter_names=parameter_names,
parameter_units=parameter_units,
- bounds=[D_bounds, slope_bounds],
+ priors=priors,
coordinate_name='box_length',
)
@@ -159,9 +161,35 @@ def linear_func(x, a, b):
# Use these as initial parameters for optimization
x0 = [D_0_init, slope_init]
- # Convert bounds to format expected by scipy
- bounds_scipy = [(b[0].value, b[1].value) for b in self.bounds]
- result = minimize(self.nll, x0, bounds=bounds_scipy, method='L-BFGS-B')
+ result = minimize(self.nll, x0, method='L-BFGS-B')
+
+ # Store results
+ self.data_group['D_0'] = result.x[0] * self.parameter_units[0]
+ self.data_group['slope'] = result.x[1] * self.parameter_units[1]
+
+ def max_aposteriori(self):
+ """Find maximum aposteriori parameters with better initial guess for YehHummer."""
+ # Use linear fit for initial parameters
+ inv_L, D_values, D_errors = self._prepare_data_for_fit()
+
+ def linear_func(x, a, b):
+ return a - b * x
+
+ popt, _ = curve_fit(
+ linear_func,
+ inv_L,
+ D_values,
+ sigma=D_errors if np.any(D_errors > 0) else None,
+ p0=[np.max(D_values), (D_values[0] - D_values[-1]) / (inv_L[0] - inv_L[-1])],
+ )
+
+ D_0_init = popt[0]
+ slope_init = popt[1]
+
+ # Use these as initial parameters for optimization
+ x0 = [D_0_init, slope_init]
+
+ result = minimize(self.nlp, x0, method='L-BFGS-B')
# Store results
self.data_group['D_0'] = result.x[0] * self.parameter_units[0]