diff --git a/.cirrus.yml b/.cirrus.yml
index c38b727..579b09c 100644
--- a/.cirrus.yml
+++ b/.cirrus.yml
@@ -54,24 +54,27 @@ build_and_test_task:
     folder: ~/.cache/pip
     fingerprint_script: echo $PYTHON_VERSION
   before_script: |
-    python -m pip install --upgrade pip
+    python -m uv pip install --upgrade pip uv
   build_script: |
-    python -m pip install '.[test]' --verbose
+    python -m uv pip install '.[test]' --verbose
   test_script: |
     python -m pytest tests/ --cov=primate --benchmark-skip  
+    coverage report -m
   coverage_script: |
     pipx run 'coverage[toml]' xml -o coverage.xml --rcfile=pyproject.toml
     pipx run 'coveralls<4' --submit coverage.xml --rcfile=pyproject.toml
     pipx run 'coveralls<4' --finish
-  uninstall_script: |
-    python -m pip uninstall primate --yes
-  wheel_script: |
-    python -m build
-  install_script: |
-    python -m pip install dist/scikit_primate*.whl
-  test_coverage_script: |
-    python -m pytest tests/ --cov=primate --benchmark-skip
-    coverage report -m
+
+
+# uninstall_script: |
+#   python -m pip uninstall primate --yes
+# wheel_script: |
+#   python -m build
+# install_script: |
+#   python -m pip install dist/scikit_primate*.whl
+# test_coverage_script: |
+#   python -m pytest tests/ --cov=primate --benchmark-skip
+#   coverage report -m
 
 # coverage_report_task:
 #   # only_if: changesInclude('.cirrus.yml', '**.{h,cpp,py,toml,build}')
diff --git a/dagger.json b/dagger.json
new file mode 100644
index 0000000..6790b5c
--- /dev/null
+++ b/dagger.json
@@ -0,0 +1,4 @@
+{
+  "name": "primate",
+  "engineVersion": "v0.14.0"
+}
diff --git a/dagger.py b/dagger.py
new file mode 100644
index 0000000..4e0562b
--- /dev/null
+++ b/dagger.py
@@ -0,0 +1,20 @@
+import sys
+
+import anyio
+
+# import dagger
+from dagger import dag, connection
+
+
+async def main(args: list[str]):
+	async with connection():
+		# build container with cowsay entrypoint
+		ctr = dag.container().from_("python:alpine").with_exec(["pip", "install", "cowsay"]).with_entrypoint(["cowsay"])
+
+		# run cowsay with requested message
+		result = await ctr.with_exec(args).stdout()
+
+	print(result)
+
+
+anyio.run(main, sys.argv[1:])
diff --git a/src/primate/diagonal.py b/src/primate/diagonal.py
index 3d515f3..8f2543f 100644
--- a/src/primate/diagonal.py
+++ b/src/primate/diagonal.py
@@ -57,17 +57,18 @@ def diag(
 	## Parameterize the random vector generation
 	rng = np.random.default_rng(seed)
 	pdf = isotropic(pdf=pdf, seed=rng)
-	estimator = MeanEstimator(kwargs.pop("record", False))
+	estimator = MeanEstimator(record=kwargs.pop("record", False))
 	converge = convergence_criterion(converge, **kwargs)
 
 	## Catch degenerate case
 	if np.prod(A.shape) == 0:
-		return 0.0 if not full else (0.0, EstimatorResult(0.0, False, "", 0, []))
+		return 0.0 if not full else (0.0, EstimatorResult())
 
 	## Commence the Monte-Carlo iterations
 	if full or callback is not None:
 		numer, denom = np.zeros(N, dtype=f_dtype), np.zeros(N, dtype=f_dtype)
-		result = EstimatorResult(numer, False, converge, 0, {})
+		result = EstimatorResult(estimator, converge)
+
 		while not converge(estimator):
 			v = pdf(size=(N, 1), seed=rng).astype(f_dtype)
 			u = (A @ v).ravel()
diff --git a/src/primate/estimators.py b/src/primate/estimators.py
index e9f3f98..78c10f3 100644
--- a/src/primate/estimators.py
+++ b/src/primate/estimators.py
@@ -1,6 +1,9 @@
+"""Estimators API."""
+
 from dataclasses import dataclass, field
+from numbers import Integral
 from operator import and_, or_
-from typing import Callable, Optional, Protocol, Sized, Union, runtime_checkable
+from typing import Callable, Iterable, Literal, Optional, Protocol, Sized, Union, runtime_checkable
 
 import numpy as np
 import scipy as sp
@@ -8,7 +11,7 @@
 from .stats import Covariance
 
 
-def arr_summary(x: Union[float, np.ndarray]):
+def arr_summary(x: Union[float, np.ndarray]) -> str:
 	if x is None:
 		return "None"
 	x = np.atleast_1d(x)
@@ -37,7 +40,7 @@ def update(self, x: Union[float, np.ndarray], **kwargs: dict): ...
 	def estimate(self) -> Union[float, np.ndarray]: ...
 
 
-class ConvergenceCriterion(Callable):
+class ConvergenceCriterion:
 	"""Generic stopping criteria for sequences."""
 
 	def __init__(self, operation: Callable):
@@ -58,40 +61,50 @@ def __call__(self, est: Estimator) -> bool:
 
 @dataclass
 class EstimatorResult:
-	estimate: Union[float, np.ndarray]
-	estimator: Estimator
+	"""Data class for representing the results of statistical estimators."""
+
+	estimator: Optional[Estimator] = None
 	criterion: Union[ConvergenceCriterion, str, None] = None
+	estimate: Union[float, np.ndarray] = 0.0
 	message: str = ""
 	nit: int = 0
 	info: dict = field(default_factory=dict)
 
+	def __iter__(self) -> Iterable:
+		return iter((self.estimator, self.criterion, self.estimate, self.message, self.nit, self.info))
+
 	def update(self, est: Estimator, converge: ConvergenceCriterion, **kwargs: dict):
-		self.estimate = est.estimate
 		self.estimator = est
+		self.criterion = converge
+		self.estimate = est.estimate
 		self.message = converge.message(est) if hasattr(converge, "message") else ""
-		self.nit = est.__len__()
-		self.info = self.info | kwargs
+		self.nit = len(est)
+		self.info |= kwargs
 
 
 class MeanEstimator(Estimator):
 	"""Sample mean estimator with stable covariance updating."""
 
-	delta: Union[float, np.ndarray]
-	cov: Optional[Covariance]
-	values: Optional[list]
+	n_samples: int = 0
+	delta: Union[float, np.ndarray] = np.atleast_1d(np.inf)
+	cov: Optional[Covariance] = None
+	values: Optional[list] = []
 
-	def __init__(self, record: bool = False) -> None:
+	def __init__(self, dim: Optional[int] = None, record: bool = False) -> None:
 		self.n_samples = 0
-		self.delta = None
-		self.cov = None
+		self.delta = np.atleast_1d(np.inf)
+		self.cov = Covariance(dim=dim) if isinstance(dim, Integral) else None
 		self.values = [] if record else None
 
 	@property
 	def mean(self) -> Optional[float]:
 		if self.cov is None:
 			return None
-		return self.cov.mu.item() if len(self.cov.mu) == 1 else self.cov.mu.ravel()
+		mu = np.atleast_1d(self.cov.mean)
+		return mu.item() if len(mu) == 1 else np.ravel(mu)
 
+	## ndim = 1, x.shape = (n,) => we have n samples of 1-d
+	## ndim = 2, x.shape = (n,m) => we have n samples of m-d
 	def update(self, x: Union[float, np.ndarray]):
 		x = np.atleast_1d(x)
 		x = x[:, None] if x.ndim == 1 else x
@@ -179,7 +192,9 @@ def message(self, est: MeanEstimator) -> bool:
 
 
 class ToleranceCriterion(ConvergenceCriterion):
-	def __init__(self, rtol: float = 0.01, atol: float = 1.49e-08, ord: Union[int, str] = 2) -> None:
+	def __init__(
+		self, rtol: float = 0.01, atol: float = 1.49e-08, ord: Union[Literal["fro", "nuc"], float, None] = 2.0
+	) -> None:
 		self.rtol = rtol
 		self.atol = atol
 		self.ord = ord
@@ -188,14 +203,15 @@ def __call__(self, est: MeanEstimator) -> bool:
 		if est.mean is None:
 			return False
 		error = np.linalg.norm(est.delta, ord=self.ord)
-		return error < self.atol or error < self.rtol * np.linalg.norm(est.mean, ord=self.ord)
+		estimate = np.atleast_1d(est.estimate)
+		return bool(error < self.atol or error < self.rtol * np.linalg.norm(estimate, ord=self.ord))
 
 	def message(self, est: MeanEstimator) -> str:
 		msg = f"Est: {arr_summary(est.estimate)}"
 		msg += f"(atol={self.atol:3f}, rtol={self.rtol:3f}, #S:{ len(est) })"
 		if est.estimate is not None:
-			error = np.linalg.norm(est.delta, ord=self.ord)
-			norm = np.linalg.norm(est.estimate, ord=self.ord)
+			error = np.linalg.norm(np.atleast_1d(est.delta), ord=self.ord)
+			norm = np.linalg.norm(np.atleast_1d(est.estimate), ord=self.ord)
 			msg += f"\nnorm(it - est, {self.ord}) = {error:.3f}, norm(est, {self.ord}) = {norm:.3f}"
 		return msg
 
@@ -240,7 +256,7 @@ def __init__(self, confidence: float = 0.95, atol: float = 0.05, rtol: float = 0
 	# 	self.vr_est = L * self.vr_pre + denom * (estimate - self.mu_pre) ** 2  # update sample variance
 	# 	self.mu_pre = self.mu_est
 	# 	self.vr_pre = self.vr_est
-	def _error(self, est: Estimator):
+	def _error(self, est: MeanEstimator):
 		if est.n_samples < 3:
 			return (np.inf, np.inf)
 		std_dev = est.cov.covariance() ** (1 / 2)
@@ -265,7 +281,7 @@ class KneeCriterion(ConvergenceCriterion):
 	def __init__(self, S: float = 1.0) -> None:
 		self.S = S
 
-	def __call__(self, est: MeanEstimator):
+	def __call__(self, est: MeanEstimator) -> bool:
 		"""Applies the kneedle algorithm to detect the knee in the sequence."""
 		if est.values is None or len(est.values) < 3:
 			return False
@@ -303,7 +319,8 @@ def message(self, est: MeanEstimator) -> str:
 }
 
 
-def convergence_criterion(criterion: Union[str, ConvergenceCriterion], **kwargs) -> ConvergenceCriterion:
+def convergence_criterion(criterion: Union[str, ConvergenceCriterion], **kwargs: dict) -> ConvergenceCriterion:
+	"""Parameterizes a convergence criterion."""
 	# assert criterion.lower() in CRITERIA.keys(), f"Invalid criterion {criterion}"
 	if isinstance(criterion, ConvergenceCriterion):
 		return criterion
diff --git a/src/primate/random.py b/src/primate/random.py
index 58578a1..bad5175 100644
--- a/src/primate/random.py
+++ b/src/primate/random.py
@@ -16,7 +16,7 @@
 def symmetric(
 	n: int,
 	dist: str = "normal",
-	pd: bool = True,
+	pd: bool = False,
 	ew: Optional[np.ndarray] = None,
 	seed: Union[int, np.random.Generator, None] = None,
 ) -> np.ndarray:
@@ -73,7 +73,9 @@ def haar(n: int, ew: Optional[np.ndarray] = None, seed: Union[int, np.random.Gen
 
 
 def isotropic(
-	size: Union[int, tuple, None] = None, pdf: str = "rademacher", seed: Union[int, np.random.Generator, None] = None
+	size: Union[int, tuple, None] = None,
+	pdf: Union[Callable, str] = "rademacher",
+	seed: Union[int, np.random.Generator, None] = None,
 ) -> Union[np.ndarray, Callable]:
 	"""Generates random vectors from a specified isotropic distribution.
 
@@ -85,6 +87,8 @@ def isotropic(
 	Returns:
 		Array of shape `size` with rows distributed according to `pdf`.
 	"""
+	if callable(pdf):
+		return pdf
 	assert pdf in _ISO_DISTRIBUTIONS.keys(), f"Invalid distribution '{pdf}' supplied."
 	pdf: str = _ISO_DISTRIBUTIONS[pdf]
 	rng = np.random.default_rng(seed)
diff --git a/src/primate/stats.py b/src/primate/stats.py
index 749fcd3..560d0dc 100644
--- a/src/primate/stats.py
+++ b/src/primate/stats.py
@@ -17,6 +17,8 @@ def __init__(self, dim: int = 1):
 
 	@property
 	def mean(self):
+		if self.n == 0:
+			return np.nan
 		return self.mu.item() if self.dim == 1 else self.mu
 
 	def update(self, X: np.ndarray) -> None:
diff --git a/src/primate/trace.py b/src/primate/trace.py
index 2b1fe3b..e640e27 100644
--- a/src/primate/trace.py
+++ b/src/primate/trace.py
@@ -113,18 +113,18 @@ def hutch(
 	## Parameterize the various quantities
 	rng = np.random.default_rng(seed)
 	pdf = isotropic(pdf=pdf, seed=rng)
-	estimator = MeanEstimator(kwargs.pop("record", False))
+	estimator = MeanEstimator(record=kwargs.pop("record", False))
 	converge = convergence_criterion(converge, **kwargs)
 	# quad_form = (lambda v: A.quad(v)) if hasattr(A, "quad") else (lambda v: (v.T @ (A @ v)).item())
 	quad_form = (lambda v: A.quad(v)) if hasattr(A, "quad") else (lambda v: np.diag(np.atleast_2d((v.T @ (A @ v)))))
 
 	## Catch degenerate case
 	if np.prod(A.shape) == 0:
-		return 0.0 if not full else (0.0, EstimatorResult(0.0, False, converge, 0, {}))
+		return 0.0 if not full else (0.0, EstimatorResult(estimator, converge))
 
 	## Commence the Monte-Carlo iterations
 	if full or callback is not None:
-		result = EstimatorResult(0.0, False, converge, 0, {})
+		result = EstimatorResult(estimator, converge)
 		callback = lambda x: x if callback is None else callback
 		while not converge(estimator):
 			v = pdf(size=(N, batch)).astype(f_dtype)
@@ -287,7 +287,7 @@ def xtrace(
 	## Commence the batched-iteration
 	estimate = np.inf
 	it = 0
-	result = EstimatorResult(0.0, False, None, 0, {})
+	result = EstimatorResult()
 	rng = np.random.default_rng(seed)
 	while (it * batch) < A.shape[1]:  # err >= (error_atol + error_rtol * abs(t)):
 		## Determine number of new sample vectors to generate
diff --git a/tests/test_estimators.py b/tests/test_estimators.py
index 0203bbb..7bf04ca 100644
--- a/tests/test_estimators.py
+++ b/tests/test_estimators.py
@@ -66,7 +66,7 @@ def test_CountCriterion():
 
 def test_ToleranceCriterion():
 	rng = np.random.default_rng(1234)
-	mu = MeanEstimator()
+	mu = MeanEstimator(15)
 	cc = ToleranceCriterion(atol=0, rtol=0.10, ord=1)
 	while not cc(mu):
 		mu.update(rng.uniform(size=(1, 15), low=-1, high=+1))
diff --git a/tests/test_operator.py b/tests/test_operator.py
index c8663f8..c888b17 100644
--- a/tests/test_operator.py
+++ b/tests/test_operator.py
@@ -69,10 +69,10 @@ def test_operator_interface():
 	assert np.allclose(A @ v, M @ v)
 
 
-def test_spectral():
+def test_spectral_positive_definite():
 	rng = np.random.default_rng(1234)
 	n = 100
-	A = symmetric(n)
+	A = symmetric(n, pd=True)
 	v = rng.uniform(size=A.shape[1], low=-1, high=1)
 	for fun in _BUILTIN_MATRIX_FUNCTIONS:
 		f = param_callable(fun)
diff --git a/tests/test_quadrature.py b/tests/test_quadrature.py
index 1c43022..198e805 100644
--- a/tests/test_quadrature.py
+++ b/tests/test_quadrature.py
@@ -6,7 +6,7 @@
 
 def test_quadrature():
 	rng = np.random.default_rng(seed=1234)
-	A = symmetric(50, seed=rng)
+	A = symmetric(50, seed=rng, pd=True)
 	quad_ests = []
 	for _ in range(100):
 		v = rng.uniform(size=A.shape[1], low=0, high=1)