ParticleGroup.cov() is implemented as:
def cov(self, *keys):
"""
Covariance matrix from any properties
"""
dats = np.array([self[key] for key in keys])
return np.cov(dats, aweights=self.weight)
On numpy ≥ 2.x, np.cov defaults to ddof=1 (i.e., bias=False) regardless of whether aweights is provided. This means the covariance is divided by N-1 introducing a systematic 1/(N-1) relative bias.
For typical use cases with large N this is negligible, but it becomes visible when verifying exact covariance reproduction (e.g., from a known input matrix). With N=200,000 particles, the diagonal elements are inflated by 5×10⁻⁶ relative.
Demonstration via Distgen:
import numpy as np
from distgen import Generator
D = Generator('gaussian_nd.in.yaml', verbose=0)
P = D.run()
coords = ['x', 'y', 'z', 'px', 'py', 'pz']
cov_input = np.array(D._input['nd_gaussian_dist']['cov_matrix'], dtype=float)
# P.cov() uses ddof=1 implicitly
diag_out = np.diag(P.cov(*coords))
diag_in = np.diag(cov_input)
print('P.cov() rel error:', np.max(np.abs((diag_out - diag_in)/diag_in)))
# → 5.0e-06 = 1/(N-1)
# Manual ddof=0 shows machine precision
data = np.array([getattr(P, c) for c in coords])
diag_ddof0 = np.diag(np.cov(data, ddof=0))
print('ddof=0 rel error:', np.max(np.abs((diag_ddof0 - diag_in)/diag_in)))
# → ~1e-15
Note: In older numpy versions (< 2.0), the documentation stated that ddof defaults to 0 when aweights or fweights are provided. This appears to have changed in numpy 2.x where bias=False (default) always sets ddof=1.
Suggestion: Consider adding a ddof parameter to ParticleGroup.cov(), or documenting the current behavior, or passing bias=True to get ddof=0 (population covariance).
ParticleGroup.cov() is implemented as:
On numpy ≥ 2.x,
np.covdefaults toddof=1(i.e., bias=False) regardless of whether aweights is provided. This means the covariance is divided by N-1 introducing a systematic1/(N-1)relative bias.For typical use cases with large N this is negligible, but it becomes visible when verifying exact covariance reproduction (e.g., from a known input matrix). With N=200,000 particles, the diagonal elements are inflated by 5×10⁻⁶ relative.
Demonstration via Distgen:
Note: In older numpy versions (< 2.0), the documentation stated that ddof defaults to 0 when aweights or fweights are provided. This appears to have changed in numpy 2.x where bias=False (default) always sets ddof=1.
Suggestion: Consider adding a ddof parameter to ParticleGroup.cov(), or documenting the current behavior, or passing bias=True to get ddof=0 (population covariance).