diff --git a/n3fit/src/n3fit/tests/test_vpinterface.py b/n3fit/src/n3fit/tests/test_vpinterface.py
index 5a90d6cbf9..26e77ce2cd 100644
--- a/n3fit/src/n3fit/tests/test_vpinterface.py
+++ b/n3fit/src/n3fit/tests/test_vpinterface.py
@@ -56,7 +56,7 @@ def test_N3PDF(members, layers):
     assert n3pdf.compute_arclength().shape == (5,)
     # Try to get a plotting grid
     res = xplotting_grid(n3pdf, 1.6, xx)
-    assert res.grid_values.shape == (members, 8, xsize)
+    assert res.grid_values.data.shape == (members, 8, xsize)
 
 
 def test_vpinterface():
@@ -68,7 +68,7 @@ def test_vpinterface():
     res_2 = xplotting_grid(fit_2, 1.6, xgrid)
     distances = distance_grids([fit_1, fit_2], [res_1, res_2], 0)
     assert len(distances) == 2
-    assert distances[0].grid_values.shape == (8, 40)
-    assert distances[1].grid_values.shape == (8, 40)
-    np.testing.assert_allclose(distances[0].grid_values, 0.0)
-    assert not np.allclose(distances[1].grid_values, 0.0)
+    assert distances[0].grid_values.data.shape == (8, 40)
+    assert distances[1].grid_values.data.shape == (8, 40)
+    np.testing.assert_allclose(distances[0].grid_values.data, 0.0)
+    assert not np.allclose(distances[1].grid_values.data, 0.0)
diff --git a/n3fit/src/n3fit/vpinterface.py b/n3fit/src/n3fit/vpinterface.py
index 28af0a0469..98ccd92bed 100644
--- a/n3fit/src/n3fit/vpinterface.py
+++ b/n3fit/src/n3fit/vpinterface.py
@@ -13,7 +13,7 @@
     >>> pdf_model = pdfNN_layer_generator(nodes=[8], activations=['linear'], seed=0, flav_info=fake_fl)
     >>> n3pdf = N3PDF(pdf_model)
     >>> res = xplotting_grid(n3pdf, 1.6, fake_x)
-    >>> res.grid_values.shape
+    >>> res.grid_values.error_members().shape
     (1, 8, 3)
 
 
diff --git a/validphys2/src/validphys/arclength.py b/validphys2/src/validphys/arclength.py
index 5ae170fc25..3484c5742e 100644
--- a/validphys2/src/validphys/arclength.py
+++ b/validphys2/src/validphys/arclength.py
@@ -50,7 +50,7 @@ def arc_lengths(
         eps = (b - a) / npoints
         ixgrid = xgrid(a, b, "linear", npoints)
         # PDFs evaluated on grid
-        xfgrid = xplotting_grid(pdf, Q, ixgrid, basis, flavours).grid_values * ixgrid[1]
+        xfgrid = xplotting_grid(pdf, Q, ixgrid, basis, flavours).grid_values.data * ixgrid[1]
         fdiff = np.diff(xfgrid) / eps  # Compute forward differences
         res += integrate.simps(1 + np.square(fdiff), ixgrid[1][1:])
     stats = pdf.stats_class(res)
@@ -127,6 +127,6 @@ def integrability_number(
     checked = check_basis(basis, flavours)
     basis, flavours = checked["basis"], checked["flavours"]
     ixgrid = xgrid(1e-9, 1e-6, "log", 3)
-    xfgrid = xplotting_grid(pdf, Q, ixgrid, basis, flavours).grid_values
+    xfgrid = xplotting_grid(pdf, Q, ixgrid, basis, flavours).grid_values.data
     res = np.sum(np.abs(xfgrid), axis=2)
     return res.squeeze()
diff --git a/validphys2/src/validphys/calcutils.py b/validphys2/src/validphys/calcutils.py
index af4e721735..2601e084d6 100644
--- a/validphys2/src/validphys/calcutils.py
+++ b/validphys2/src/validphys/calcutils.py
@@ -70,7 +70,7 @@ def all_chi2(results):
     """Return the chi² for all elements in the result. Note that the
     interpretation of the result will depend on the PDF error type"""
     data_result, th_result = results
-    diffs = th_result._rawdata - data_result.central_value[:,np.newaxis]
+    diffs = th_result.rawdata - data_result.central_value[:,np.newaxis]
     return calc_chi2(sqrtcov=data_result.sqrtcovmat, diffs=diffs)
 
 def central_chi2(results):
@@ -85,7 +85,7 @@ def all_chi2_theory(results, totcov):
     """Like all_chi2 but here the chi² are calculated using a covariance matrix
     that is the sum of the experimental covmat and the theory covmat."""
     data_result, th_result = results
-    diffs = th_result._rawdata - data_result.central_value[:,np.newaxis]
+    diffs = th_result.rawdata - data_result.central_value[:,np.newaxis]
     total_covmat = np.array(totcov)
     return calc_chi2(sqrtcov=la.cholesky(total_covmat, lower=True), diffs=diffs)
 
diff --git a/validphys2/src/validphys/closuretest/closure_results.py b/validphys2/src/validphys/closuretest/closure_results.py
index 9f30cf15ee..857720fe4a 100644
--- a/validphys2/src/validphys/closuretest/closure_results.py
+++ b/validphys2/src/validphys/closuretest/closure_results.py
@@ -122,7 +122,7 @@ def bootstrap_bias_experiment(
     """
     dt_ct, th_ct = dataset_inputs_results
     ((_, th_ul),) = underlying_dataset_inputs_results
-    th_ct_boot_cv = bootstrap_values(th_ct._rawdata, bootstrap_samples)
+    th_ct_boot_cv = bootstrap_values(th_ct.rawdata, bootstrap_samples)
     boot_diffs = th_ct_boot_cv - th_ul.central_value[:, np.newaxis]
     boot_bias = calc_chi2(dt_ct.sqrtcovmat, boot_diffs) / len(dt_ct)
     return boot_bias
@@ -191,7 +191,7 @@ def variance_dataset(results, fit, use_fitcommondata):
 
     """
     dt, th = results
-    diff = th.central_value[:, np.newaxis] - th._rawdata
+    diff = th.central_value[:, np.newaxis] - th.rawdata
     var_unnorm = calc_chi2(dt.sqrtcovmat, diff).mean()
     return VarianceData(var_unnorm, len(th))
 
@@ -210,7 +210,7 @@ def bootstrap_variance_experiment(dataset_inputs_results, bootstrap_samples=500)
     normalised to the number of data in the experiment.
     """
     dt_ct, th_ct = dataset_inputs_results
-    diff = th_ct.central_value[:, np.newaxis] - th_ct._rawdata
+    diff = th_ct.central_value[:, np.newaxis] - th_ct.rawdata
     var_unnorm_boot = bootstrap_values(
         diff,
         bootstrap_samples,
diff --git a/validphys2/src/validphys/closuretest/multiclosure.py b/validphys2/src/validphys/closuretest/multiclosure.py
index c784c8996a..c48388b238 100644
--- a/validphys2/src/validphys/closuretest/multiclosure.py
+++ b/validphys2/src/validphys/closuretest/multiclosure.py
@@ -130,7 +130,7 @@ def fits_dataset_bias_variance(
     """
     closures_th, law_th, _, sqrtcov = internal_multiclosure_dataset_loader
     # The dimentions here are (fit, data point, replica)
-    reps = np.asarray([th._rawdata[:, :_internal_max_reps] for th in closures_th])
+    reps = np.asarray([th.rawdata[:, :_internal_max_reps] for th in closures_th])
     # take mean across replicas - since we might have changed no. of reps
     centrals = reps.mean(axis=2)
     # place bins on first axis
@@ -226,7 +226,7 @@ def dataset_xi(internal_multiclosure_dataset_loader):
 
     """
     closures_th, law_th, covmat, _ = internal_multiclosure_dataset_loader
-    replicas = np.asarray([th._rawdata for th in closures_th])
+    replicas = np.asarray([th.rawdata for th in closures_th])
     centrals = np.mean(replicas, axis=-1)
     underlying = law_th.central_value
 
@@ -297,7 +297,7 @@ class BootstrappedTheoryResult:
     """
 
     def __init__(self, data):
-        self._rawdata = data
+        self.rawdata = data
         self.central_value = data.mean(axis=1)
 
 
@@ -341,7 +341,7 @@ def _bootstrap_multiclosure_fits(
     # construct proxy fits theory predictions
     for fit_th in fit_boot_th:
         rep_boot_index = rng.choice(n_rep_max, size=n_rep, replace=use_repeats)
-        boot_ths.append(BootstrappedTheoryResult(fit_th._rawdata[:, rep_boot_index]))
+        boot_ths.append(BootstrappedTheoryResult(fit_th.rawdata[:, rep_boot_index]))
     return (boot_ths, *input_tuple)
 
 
@@ -741,7 +741,7 @@ def dataset_fits_bias_replicas_variance_samples(
     """
     closures_th, law_th, _, sqrtcov = internal_multiclosure_dataset_loader
     # The dimentions here are (fit, data point, replica)
-    reps = np.asarray([th._rawdata[:, :_internal_max_reps] for th in closures_th])
+    reps = np.asarray([th.rawdata[:, :_internal_max_reps] for th in closures_th])
     # take mean across replicas - since we might have changed no. of reps
     centrals = reps.mean(axis=2)
     # place bins on first axis
diff --git a/validphys2/src/validphys/closuretest/multiclosure_pdf.py b/validphys2/src/validphys/closuretest/multiclosure_pdf.py
index 7881c1e557..a8d0f5f3f8 100644
--- a/validphys2/src/validphys/closuretest/multiclosure_pdf.py
+++ b/validphys2/src/validphys/closuretest/multiclosure_pdf.py
@@ -90,7 +90,7 @@ def xi_grid_values(xi_pdfgrids):
     glu_sin_grid, nonsin_grid = xi_pdfgrids
     # grid values have shape: replica, flavour, x
     # concatenate along flavour
-    return np.concatenate((glu_sin_grid.grid_values, nonsin_grid.grid_values), axis=1)
+    return np.concatenate((glu_sin_grid.grid_values.data, nonsin_grid.grid_values.data), axis=1)
 
 
 def underlying_xi_grid_values(
diff --git a/validphys2/src/validphys/closuretest/multiclosure_pseudodata.py b/validphys2/src/validphys/closuretest/multiclosure_pseudodata.py
index ff32d86fd1..22b2c78a1a 100644
--- a/validphys2/src/validphys/closuretest/multiclosure_pseudodata.py
+++ b/validphys2/src/validphys/closuretest/multiclosure_pseudodata.py
@@ -126,7 +126,7 @@ def experiments_closure_pseudodata_estimators_table(
         fits_erep_delta_chi2 = []
         for i_fit, fit_th in enumerate(closures_th):
             # some of these could be done outside of loop, but easier to do here.
-            th_replicas = fit_th._rawdata
+            th_replicas = fit_th.rawdata
             th_central = np.mean(th_replicas, axis=-1)
             dt_replicas = fits_reps_pseudo[i_fit].xs(exp_name, axis=0, level=0).to_numpy()
             dt_central = fits_cv.xs(exp_name, axis=0, level=0).iloc[:, i_fit].to_numpy()
diff --git a/validphys2/src/validphys/correlations.py b/validphys2/src/validphys/correlations.py
index 0691a71f96..2e6132c3cf 100644
--- a/validphys2/src/validphys/correlations.py
+++ b/validphys2/src/validphys/correlations.py
@@ -14,12 +14,12 @@
 
 #This would be a good candidate to be optimized to calculate everything in one
 #pass over x,
-def _basic_obs_pdf_correlation(pdfarr, obsarr):
+def _basic_obs_pdf_correlation(pdf_stats, obs_stats):
     """Calculate the correlation between pdfs and observables.
-    The expected format is:
+    The expected format is are Stats instances of:
 
-    obsarr: (nbins x nreplicas), as returned from thresult.
-    pdfarr: (nreplicas x nf x nx), as returned from xplotting_grid.grid_values
+    obs_stats: (nbins x nreplicas), as returned from thresult.
+    pdf_stats: (nreplicas x nf x nx), as returned from xplotting_grid.grid_values
 
     The returned array contains the PDF correlation between
     the value of the obsevable and the PDF at the corresponding point in (fl,x)
@@ -27,11 +27,8 @@ def _basic_obs_pdf_correlation(pdfarr, obsarr):
     (nbins x nf x nx), compatible with grid_values, upon
     changing replicas->bins.
     """
-
-    #Remove mean
-    #TODO: This should be done at the Result level
-    x = pdfarr  - np.mean(pdfarr, axis=0)
-    y = obsarr.T - np.mean(obsarr, axis=1)
+    x = pdf_stats.data - pdf_stats.central_value()
+    y = obs_stats.data - obs_stats.central_value()
 
     #We want to compute:
     #sum(x*y)/(norm(x)*norm(y))
@@ -46,28 +43,29 @@ def _basic_obs_pdf_correlation(pdfarr, obsarr):
 
     return num/den
 
-def _basic_obs_obs_correlation(obsarr1, obsarr2):
+def _basic_obs_obs_correlation(obs1, obs2):
     """Calculate the correlation between two observables. The expected format is
-    obsarr1: (nbins1, nreplicas)
-    obsarr2: (nbins2, nreplicas)
+    Stats instances of:
+
+    obs1: (nreplicas, nbins1)
+    obs2: (nreplicas, nbins2)
 
     The result is (nbins1 , nbins2), a mareix containing the correlation
     coefficients between the two sets.
     """
-    #TODO: Do this at Result level taking into account error type
-    x = (obsarr1.T - np.mean(obsarr1, axis=1)).T
-    y = (obsarr2.T - np.mean(obsarr2, axis=1))
+    x = (obs1.data - obs1.central_value()).T
+    y = (obs2.data - obs2.central_value())
 
     return x@y/np.outer(la.norm(x,axis=1),la.norm(y,axis=0))
 
-#TODO: Implement for other error types. Do not use the _rawdata.
+#TODO: Implement for other error types. Do not use the rawdata.
 @check_pdf_is_montecarlo
 def obs_pdf_correlations(pdf, results, xplotting_grid):
     """Return the correlations between each point in a dataset and the PDF
     values on a grid of (x,f) points in a format similar to `xplotting_grid`."""
     _ , th = results
-    corrs = _basic_obs_pdf_correlation(xplotting_grid.grid_values, th._rawdata)
-    return xplotting_grid._replace(grid_values=corrs)
+    corrs = pdf.stats_class(_basic_obs_pdf_correlation(xplotting_grid.grid_values, th.stats))
+    return xplotting_grid.copy_grid(grid_values=corrs)
 
 
 corrpair_results = collect('results', ['corrpair'])
@@ -77,4 +75,4 @@ def obs_pdf_correlations(pdf, results, xplotting_grid):
 def obs_obs_correlations(pdf, corrpair_results):
     """Return the theoretical correlation matrix between a pair of observables."""
     (_,th1), (_,th2) = corrpair_results
-    return _basic_obs_obs_correlation(th1._rawdata, th2._rawdata)
+    return _basic_obs_obs_correlation(th1.stats, th2.stats)
diff --git a/validphys2/src/validphys/covmats.py b/validphys2/src/validphys/covmats.py
index bd09c15d5b..d4d389ea6b 100644
--- a/validphys2/src/validphys/covmats.py
+++ b/validphys2/src/validphys/covmats.py
@@ -531,7 +531,7 @@ def pdferr_plus_covmat(dataset, pdf, covmat_t0_considered):
     True
     """
     th = ThPredictionsResult.from_convolution(pdf, dataset)
-    pdf_cov = np.cov(th._rawdata, rowvar=True)
+    pdf_cov = np.cov(th.rawdata, rowvar=True)
     return pdf_cov + covmat_t0_considered
 
 
diff --git a/validphys2/src/validphys/dataplots.py b/validphys2/src/validphys/dataplots.py
index bbb1279f50..b9bfd1def3 100644
--- a/validphys2/src/validphys/dataplots.py
+++ b/validphys2/src/validphys/dataplots.py
@@ -871,7 +871,7 @@ def plot_smpdf(pdf, dataset, obs_pdf_correlations, mark_threshold:float=0.9):
 
     basis = obs_pdf_correlations.basis
 
-    fullgrid = obs_pdf_correlations.grid_values
+    fullgrid = obs_pdf_correlations.grid_values.data
 
     fls = obs_pdf_correlations.flavours
     x = obs_pdf_correlations.xgrid
diff --git a/validphys2/src/validphys/deltachi2.py b/validphys2/src/validphys/deltachi2.py
index 65f0a9f8a8..bba69a87e5 100644
--- a/validphys2/src/validphys/deltachi2.py
+++ b/validphys2/src/validphys/deltachi2.py
@@ -134,11 +134,6 @@ def plot_delta_chi2_hessian_distribution(delta_chi2_hessian, pdf, total_chi2_dat
     return fig
 
 
-XPlottingGrid = namedtuple(
-    "XPlottingGrid", ("Q", "basis", "flavours", "xgrid", "grid_values", "scale")
-)
-
-
 def pos_neg_xplotting_grids(delta_chi2_hessian, xplotting_grid):
     """
     Generates xplotting_grids correspodning to positive and negative delta chi2s.
@@ -149,11 +144,11 @@ def pos_neg_xplotting_grids(delta_chi2_hessian, xplotting_grid):
     pos_mask = np.append(True, positive_eigenvalue_mask)
     neg_mask = np.append(True, ~positive_eigenvalue_mask)
 
-    pos_grid = xplotting_grid.grid_values[pos_mask]
-    neg_grid = xplotting_grid.grid_values[neg_mask]
+    pos_grid = xplotting_grid.grid_values.data[pos_mask]
+    neg_grid = xplotting_grid.grid_values.data[neg_mask]
 
-    pos_xplotting_grid = xplotting_grid._replace(grid_values=pos_grid)
-    neg_xplotting_grid = xplotting_grid._replace(grid_values=neg_grid)
+    pos_xplotting_grid = xplotting_grid.copy_grid(grid_values=pos_grid)
+    neg_xplotting_grid = xplotting_grid.copy_grid(grid_values=neg_grid)
 
     return [xplotting_grid, pos_xplotting_grid, neg_xplotting_grid]
 
@@ -220,7 +215,7 @@ def draw(self, pdf, grid, flstate):
         # Basically stats is an object which says what is the type class of the replicas:
         # MCStats, HessianStats, SymmHessianStats. In this way validphys use the right methods
         # to compute statistical values.
-        stats = pdf.stats_class(grid.grid_values[:, flindex, :])
+        stats = pdf.stats_class(grid.grid_values.data[:, flindex, :])
 
         # Ignore spurious normalization warnings
         with warnings.catch_warnings():
diff --git a/validphys2/src/validphys/eff_exponents.py b/validphys2/src/validphys/eff_exponents.py
index a7b07c2ecc..dce6a49d09 100644
--- a/validphys2/src/validphys/eff_exponents.py
+++ b/validphys2/src/validphys/eff_exponents.py
@@ -66,14 +66,14 @@ def alpha_eff(pdf: PDF, *,
 
     pdfGrid = pdfgrids.xplotting_grid(
         pdf, Q, xgrid=xGrid, basis=basis, flavours=flavours)
-    pdfGrid_values = pdfGrid.grid_values
+    pdfGrid_values = pdfGrid.grid_values.data
     # NOTE: without this I get "setting an array element with a sequence"
     xGrid = pdfGrid.xgrid
     with warnings.catch_warnings():
         warnings.simplefilter('ignore', RuntimeWarning)
         alphaGrid_values = -np.log(abs(pdfGrid_values/xGrid))/np.log(xGrid)
         alphaGrid_values[alphaGrid_values == - np.inf] = np.nan  # when PDF_i =0
-    alphaGrid = pdfGrid._replace(grid_values=alphaGrid_values)
+    alphaGrid = pdfGrid.copy_grid(grid_values=pdf.stats_class(alphaGrid_values))
     return alphaGrid
 
 @check_positive('Q')
@@ -110,14 +110,14 @@ def beta_eff(pdf, *,
 
     pdfGrid = pdfgrids.xplotting_grid(
         pdf, Q, xgrid=xGrid, basis=basis, flavours=flavours)
-    pdfGrid_values = pdfGrid.grid_values
+    pdfGrid_values = pdfGrid.grid_values.data
     # NOTE: without this I get "setting an array element with a sequence"
     xGrid = pdfGrid.xgrid
     with warnings.catch_warnings():
         warnings.simplefilter('ignore', RuntimeWarning)
         betaGrid_values = np.log(abs(pdfGrid_values/xGrid))/np.log(1-xGrid)
         betaGrid_values[betaGrid_values == -np.inf] = np.nan  # when PDF_i =0
-    betaGrid = pdfGrid._replace(grid_values=betaGrid_values)
+    betaGrid = pdfGrid.copy_grid(grid_values=pdf.stats_class(betaGrid_values))
 
     return betaGrid  # .grid_values
 
@@ -348,17 +348,14 @@ def next_effective_exponents_table(
 
     eff_exp_data = []
 
-    alphagrid = alpha_effs.grid_values
-    betagrid = beta_effs.grid_values
-
-    alphastats = pdf.stats_class(alphagrid)
-    betastats = pdf.stats_class(betagrid)
+    alphastats = alpha_effs.grid_values
+    betastats = beta_effs.grid_values
 
     with warnings.catch_warnings():
         warnings.simplefilter('ignore', RuntimeWarning)
 
-        alpha_cv = np.nanmean(alphagrid, axis=0)
-        beta_cv = np.nanmean(betagrid, axis=0)
+        alpha_cv = np.nanmean(alphastats.error_members(), axis=0)
+        beta_cv = np.nanmean(betastats.error_members(), axis=0)
         #tuple of low and high values repectively
         alpha68 = alphastats.errorbar68()
         beta68 = betastats.errorbar68()
diff --git a/validphys2/src/validphys/mc2hessian.py b/validphys2/src/validphys/mc2hessian.py
index 4e045ced64..2f9ee758fc 100644
--- a/validphys2/src/validphys/mc2hessian.py
+++ b/validphys2/src/validphys/mc2hessian.py
@@ -17,6 +17,8 @@
 from validphys.lhio import hessian_from_lincomb
 from validphys.pdfgrids import xplotting_grid
 
+from validphys.checks import check_pdf_is_montecarlo
+
 log = logging.getLogger(__name__)
 
 
@@ -58,6 +60,7 @@ def mc2hessian_xgrid(
     )
 
 
+@check_pdf_is_montecarlo
 def mc2hessian(
     pdf, Q, Neig: int, mc2hessian_xgrid, output_path, gridname, installgrid: bool = False
 ):
@@ -109,8 +112,8 @@ def _create_mc2hessian(pdf, Q, xgrid, Neig, output_path, name=None):
 def _get_X(pdf, Q, xgrid, reshape=False):
     pdf_grid = xplotting_grid(pdf, Q, xgrid=xgrid)
     pdf_grid_values = pdf_grid.grid_values
-    replicas = pdf_grid_values
-    mean = pdf_grid_values.mean(axis=0)
+    replicas = pdf_grid_values.error_members()
+    mean = pdf_grid_values.central_value()
     Xt = replicas - mean
     if reshape:
         Xt = Xt.reshape(Xt.shape[0], Xt.shape[1] * Xt.shape[2])
diff --git a/validphys2/src/validphys/pdfgrids.py b/validphys2/src/validphys/pdfgrids.py
index 0f9061c943..ae1eabd976 100644
--- a/validphys2/src/validphys/pdfgrids.py
+++ b/validphys2/src/validphys/pdfgrids.py
@@ -3,7 +3,9 @@
 to facilitate plotting and analysis.
 """
 from collections import namedtuple
+import dataclasses
 import numbers
+import logging
 
 import numpy as np
 import scipy.integrate as integrate
@@ -11,11 +13,13 @@
 from reportengine import collect
 from reportengine.checks import make_argcheck, CheckError, check_positive, check
 
-from validphys.core import PDF
+from validphys.core import PDF, Stats
 from validphys.gridvalues import (evaluate_luminosity)
 from validphys.pdfbases import (Basis, check_basis)
 from validphys.checks import check_pdf_normalize_to, check_xlimits
 
+log = logging.getLogger(__name__)
+
 @make_argcheck
 def _check_scale(scale):
     scales = ('linear', 'log')
@@ -37,9 +41,26 @@ def xgrid(xmin:numbers.Real=1e-5, xmax:numbers.Real=1,
     return (scale, arr)
 
 
+@dataclasses.dataclass
+class XPlottingGrid:
+    """DataClass holding the value of the PDF at the specified
+    values of x, Q and flavour.
+    The `grid_values` attribute corresponds to a `Stats` instance
+    in order to compute statistical estimators in a sensible manner.
+    """
+    Q: float
+    basis: (str, Basis)
+    flavours: (list, tuple, type(None))
+    xgrid: np.ndarray
+    grid_values: Stats
+    scale: str
 
-XPlottingGrid = namedtuple('XPlottingGrid', ('Q', 'basis', 'flavours', 'xgrid',
-                                             'grid_values', 'scale'))
+    def copy_grid(self, grid_values=None):
+        """Create a copy of the grid with potentially a different set of values"""
+        if not isinstance(grid_values, Stats):
+            log.warning("XPlottingGrid being called with a numpy grid, should be using Stats instead!")
+            grid_values = self.grid_values.__class__(grid_values)
+        return dataclasses.replace(self, grid_values=grid_values)
 
 
 @make_argcheck(check_basis)
@@ -71,8 +92,7 @@ def xplotting_grid(pdf:PDF, Q:(float,int), xgrid=None, basis:(str, Basis)='flavo
         raise TypeError(f"Invalid xgrid {xgrid!r}")
     gv = basis.grid_values(pdf, flavours, xgrid, Q)
     #Eliminante Q axis
-    #TODO: wrap this in pdf.stats_class?
-    gv = gv.reshape(gv.shape[:-1])
+    gv = pdf.stats_class(gv.reshape(gv.shape[:-1]))
 
     res = XPlottingGrid(Q, basis, flavours, xgrid, gv, scale)
     return res
@@ -233,27 +253,30 @@ def distance_grids(pdfs, xplotting_grids, normalize_to:(int,str,type(None))=None
     set is computed. At least one grid will be identical to zero.
     """
 
-    gr2 = xplotting_grids[normalize_to]
-    cv2 = pdfs[normalize_to].stats_class(gr2.grid_values).central_value()
-    sg2 = pdfs[normalize_to].stats_class(gr2.grid_values).std_error()
-    N2 = gr2.grid_values.shape[0]
+    gr2_stats = xplotting_grids[normalize_to].grid_values
+    cv2 = gr2_stats.central_value()
+    sg2 = gr2_stats.std_error()
+    N2 = pdfs[normalize_to].get_members()
 
-    newgrids = list()
+    newgrids = []
     for grid, pdf in zip(xplotting_grids, pdfs):
 
         if pdf == pdfs[normalize_to]:
-            newgrid = grid._replace(grid_values=np.zeros(shape=(grid.grid_values.shape[1], grid.grid_values.shape[2])))
+            # Zero the PDF we are normalizing against
+            pdf_zero = pdf.stats_class(np.zeros_like(gr2_stats.data[0]))
+            newgrid = grid.copy_grid(grid_values=pdf_zero)
             newgrids.append(newgrid)
             continue
 
-        cv1 = pdf.stats_class(grid.grid_values).central_value()
-        sg1 = pdf.stats_class(grid.grid_values).std_error()
-        N1 = grid.grid_values.shape[0]
+        g_stats = grid.grid_values
+        cv1 = g_stats.central_value()
+        sg1 = g_stats.std_error()
+        N1 = pdf.get_members()
 
         # the distance definition
-        distance = np.sqrt((cv1-cv2)**2/(sg1**2/N1+sg2**2/N2))
+        distance = pdf.stats_class(np.sqrt((cv1-cv2)**2/(sg1**2/N1+sg2**2/N2)))
 
-        newgrid = grid._replace(grid_values=distance)
+        newgrid = grid.copy_grid(grid_values=distance)
         newgrids.append(newgrid)
 
     return newgrids
@@ -271,29 +294,32 @@ def variance_distance_grids(pdfs, xplotting_grids, normalize_to:(int,str,type(No
     set is computed. At least one grid will be identical to zero.
     """
 
-    gr2 = xplotting_grids[normalize_to]
-    sg2 = pdfs[normalize_to].stats_class(gr2.grid_values).std_error()
-    mo2 = pdfs[normalize_to].stats_class(gr2.grid_values).moment(4)
-    N2 = gr2.grid_values.shape[0]
+    gr2_stats = xplotting_grids[normalize_to].grid_values
+    sg2 = gr2_stats.std_error()
+    mo2 = gr2_stats.moment(4)
+    N2 = pdfs[normalize_to].get_members()
     s2 = (mo2-(N2-3)/(N2-1)*sg2**4)/N2
 
-    newgrids = list()
+    newgrids = []
     for grid, pdf in zip(xplotting_grids, pdfs):
 
         if pdf == pdfs[normalize_to]:
-            newgrid = grid._replace(grid_values=np.zeros(shape=(grid.grid_values.shape[1], grid.grid_values.shape[2])))
+            # Zero the PDF we are normalizing against
+            pdf_zero = pdf.stats_class(np.zeros_like(gr2_stats.data[0]))
+            newgrid = grid.copy_grid(grid_values=pdf_zero)
             newgrids.append(newgrid)
             continue
 
-        sg1 = pdf.stats_class(grid.grid_values).std_error()
-        mo1 = pdf.stats_class(grid.grid_values).moment(4)
-        N1 = grid.grid_values.shape[0]
+        g_stats = grid.grid_values
+        sg1 = g_stats.std_error()
+        mo1 = g_stats.moment(4)
+        N1 = pdf.get_members()
         s1 = (mo1-(N1-3)/(N1-1)*sg1**4)/N1
 
         # the distance definition
-        variance_distance = np.sqrt((sg1**2-sg2**2)**2/(s1+s2))
+        variance_distance = pdf.stats_class(np.sqrt((sg1**2-sg2**2)**2/(s1+s2)))
 
-        newgrid = grid._replace(grid_values=variance_distance)
+        newgrid = grid.copy_grid(grid_values=variance_distance)
         newgrids.append(newgrid)
 
     return newgrids
diff --git a/validphys2/src/validphys/pdfplots.py b/validphys2/src/validphys/pdfplots.py
index 74397e2b83..3b668dba45 100644
--- a/validphys2/src/validphys/pdfplots.py
+++ b/validphys2/src/validphys/pdfplots.py
@@ -61,8 +61,7 @@ def normalize(self):
         if normalize_to is not None:
             normalize_pdf = self.normalize_pdf
             normalize_grid = self._xplotting_grids[normalize_to]
-            normvals = normalize_pdf.stats_class(
-                            normalize_grid.grid_values).central_value()
+            normvals = normalize_grid.grid_values.central_value()
 
             #Handle division by zero more quietly
             def fp_error(tp, flag):
@@ -75,12 +74,9 @@ def fp_error(tp, flag):
             newgrids = []
             with np.errstate(all='call'):
                 np.seterrcall(fp_error)
-                for grid in self._xplotting_grids:
-                    newvalues = grid.grid_values/normvals
-                    #newgrid is like the old grid but with updated values
-                    newgrid = type(grid)(**{**grid._asdict(),
-                                             'grid_values':newvalues})
-                    newgrids.append(newgrid)
+                for pdf, grid in zip(self.pdfs, self._xplotting_grids):
+                    newvalues = pdf.stats_class(grid.grid_values.data/normvals)
+                    newgrids.append(grid.copy_grid(grid_values=newvalues))
 
             return newgrids
         return self._xplotting_grids
@@ -190,7 +186,7 @@ def draw(self, pdf, grid, flstate):
         ax = flstate.ax
         next_prop = next(ax._get_lines.prop_cycler)
         color = next_prop['color']
-        gv = grid.grid_values[:,flstate.flindex,:]
+        gv = grid.grid_values.data[:,flstate.flindex,:]
         ax.plot(grid.xgrid, gv.T, alpha=0.2, linewidth=0.5,
                 color=color, zorder=1)
         stats = pdf.stats_class(gv)
@@ -227,7 +223,7 @@ def get_ylabel(self, parton_name):
     def draw(self, pdf, grid, flstate):
         ax = flstate.ax
         flindex = flstate.flindex
-        gv = grid.grid_values[:,flindex,:]
+        gv = grid.grid_values.data[:,flindex,:]
         stats = pdf.stats_class(gv)
 
         res = stats.std_error()
@@ -334,7 +330,7 @@ def draw(self, pdf, grid, flstate):
         next_prop = next(pcycler)
         color = next_prop['color']
 
-        gv = grid.grid_values[flindex,:]
+        gv = grid.grid_values.data[flindex,:]
         ax.plot(grid.xgrid, gv, color=color, label='$%s$' % flstate.parton_name)
 
         return gv
@@ -414,7 +410,7 @@ def draw(self, pdf, grid, flstate):
         hatchit = flstate.hatchit
         labels = flstate.labels
         handles = flstate.handles
-        stats = pdf.stats_class(grid.grid_values[:,flindex,:])
+        stats = pdf.stats_class(grid.grid_values.data[:,flindex,:])
         pcycler = ax._get_lines.prop_cycler
         #This is ugly but can't think of anything better
 
diff --git a/validphys2/src/validphys/replica_selector.py b/validphys2/src/validphys/replica_selector.py
index 200c2df36b..08a82bdf1c 100644
--- a/validphys2/src/validphys/replica_selector.py
+++ b/validphys2/src/validphys/replica_selector.py
@@ -169,7 +169,7 @@ def not_outlier_mask(
     """Return a boolean mask with the replicas that are never outside of the given
     percentile in the constrained region"""
     lim = unconstrained_region_index
-    gv = xplotting_grid.grid_values[:, :, lim:]
+    gv = xplotting_grid.grid_values.data[:, :, lim:]
     delta = nsigma * np.std(gv, axis=0)
     mean = np.mean(gv, axis=0)
     return ((gv >= mean - delta) & (gv <= mean+delta)).all(axis=2).all(axis=1)
@@ -309,7 +309,7 @@ def draw(self, pdf, grid, flstate):
             #PDFs, and instead would do something different
             next_prop = next(ax._get_lines.prop_cycler)
             color = next_prop['color']
-            gv = grid.grid_values[filt, flstate.flindex, :]
+            gv = grid.grid_values.data[filt, flstate.flindex, :]
             ax.plot(
                 grid.xgrid,
                 gv.T,
diff --git a/validphys2/src/validphys/results.py b/validphys2/src/validphys/results.py
index d989da6caa..87f95450cd 100644
--- a/validphys2/src/validphys/results.py
+++ b/validphys2/src/validphys/results.py
@@ -26,7 +26,7 @@
     check_two_dataspecs,
 )
 
-from validphys.core import DataSetSpec, PDF, DataGroupSpec
+from validphys.core import DataSetSpec, PDF, DataGroupSpec, Stats
 from validphys.calcutils import (
     all_chi2,
     central_chi2,
@@ -36,7 +36,6 @@
 )
 from validphys.convolution import (
     predictions,
-    central_predictions,
     PredictionsRequireCutsError,
 )
 
@@ -48,31 +47,15 @@ class Result:
     pass
 
 
-# TODO: Eventually,only one of (NNPDFDataResult, StatsResult) should survive
-class NNPDFDataResult(Result):
-    """A result fills its values from a pandas dataframe
-    For legacy (libNNPDF) compatibility, falls back to libNNPDF attributes"""
-
-    def __init__(self, dataobj=None, central_value=None):
-        # This class is used by both validphys and libNNPDF objects
-        # when central_value is not explictly passed, fallback to
-        # libNNPDF object .get_cv()
-        if central_value is None:
-            central_value = dataobj.get_cv()
-        self._central_value = np.array(central_value).reshape(-1)
-
-    @property
-    def central_value(self):
-        return self._central_value
-
-    def __len__(self):
-        return len(self.central_value)
-
-
 class StatsResult(Result):
     def __init__(self, stats):
         self.stats = stats
 
+    @property
+    def rawdata(self):
+        """Returns the raw data with shape (Npoints, Npdf)"""
+        return self.stats.data.T
+
     @property
     def central_value(self):
         return self.stats.central_value()
@@ -81,17 +64,28 @@ def central_value(self):
     def std_error(self):
         return self.stats.std_error()
 
+    def __len__(self):
+        """Returns the number of data points in the result"""
+        return self.rawdata.shape[0]
+
 
-class DataResult(NNPDFDataResult):
-    def __init__(self, dataobj, covmat, sqrtcovmat, central_value=None):
-        super().__init__(dataobj, central_value=central_value)
+class DataResult(StatsResult):
+    """Holds the relevant information from a given dataset"""
+    def __init__(self, dataobj, covmat, sqrtcovmat):
+        self._central_value = dataobj.get_cv()
+        stats = Stats(self._central_value)
         self._covmat = covmat
         self._sqrtcovmat = sqrtcovmat
+        super().__init__(stats)
 
     @property
     def label(self):
         return "Data"
 
+    @property
+    def central_value(self):
+        return self._central_value
+
     @property
     def std_error(self):
         return np.sqrt(np.diag(self.covmat))
@@ -106,27 +100,14 @@ def sqrtcovmat(self):
         return self._sqrtcovmat
 
 
-class ThPredictionsResult(NNPDFDataResult):
-    """Class holding theory prediction
-    For legacy purposes it still accepts libNNPDF datatypes, but prefers python-pure stuff
+class ThPredictionsResult(StatsResult):
+    """Class holding theory prediction, inherits from StatsResult
     """
-    def __init__(self, dataobj, stats_class, label=None, central_value=None):
+    def __init__(self, dataobj, stats_class, label=None):
         self.stats_class = stats_class
         self.label = label
-        # Ducktype the input into numpy arrays
-        try:
-            self._rawdata = dataobj.to_numpy()
-            # If the numpy conversion worked then we don't have a libNNPDF in our hands
-            stats = stats_class(self._rawdata.T)
-            self._std_error = stats.std_error()
-        except AttributeError:
-            self._std_error = dataobj.get_error()
-            self._rawdata = dataobj.get_data()
-        super().__init__(dataobj, central_value=central_value)
-
-    @property
-    def std_error(self):
-        return self._std_error
+        statsobj = stats_class(dataobj.T)
+        super().__init__(statsobj)
 
     @staticmethod
     def make_label(pdf, dataset):
@@ -152,20 +133,14 @@ def from_convolution(cls, pdf, dataset):
             datasets = (dataset,)
 
         try:
-            all_preds = []
-            all_centrals = []
-            for d in datasets:
-                all_preds.append(predictions(d, pdf))
-                all_centrals.append(central_predictions(d, pdf))
+            th_predictions = pd.concat([predictions(d, pdf) for d in datasets])
         except PredictionsRequireCutsError as e:
             raise PredictionsRequireCutsError("Predictions from FKTables always require cuts, "
                     "if you want to use the fktable intrinsic cuts set `use_cuts: 'internal'`") from e
-        th_predictions = pd.concat(all_preds)
-        central_values = pd.concat(all_centrals)
 
         label = cls.make_label(pdf, dataset)
 
-        return cls(th_predictions, pdf.stats_class, label, central_value=central_values)
+        return cls(th_predictions, pdf.stats_class, label)
 
 
 class PositivityResult(StatsResult):
@@ -175,10 +150,6 @@ def from_convolution(cls, pdf, posset):
         stats = pdf.stats_class(data.T)
         return cls(stats)
 
-    @property
-    def rawdata(self):
-        return self.stats.data
-
 
 # TODO: finish deprecating all dependencies on this index largely in theorycovmat module
 groups_data = collect("data", ("group_dataset_inputs_by_metadata",))
@@ -257,7 +228,7 @@ def group_result_table_no_table(groups_results, groups_index):
         ):
             replicas = (
                 ("rep_%05d" % (i + 1), th_rep)
-                for i, th_rep in enumerate(th._rawdata[index, :])
+                for i, th_rep in enumerate(th.rawdata[index, :])
             )
 
             result_records.append(
@@ -612,7 +583,7 @@ def dataset_inputs_bootstrap_phi_data(dataset_inputs_results, bootstrap_samples=
     For more information on how phi is calculated see `phi_data`
     """
     dt, th = dataset_inputs_results
-    diff = np.array(th._rawdata - dt.central_value[:, np.newaxis])
+    diff = np.array(th.rawdata - dt.central_value[:, np.newaxis])
     phi_resample = bootstrap_values(
         diff,
         bootstrap_samples,
@@ -632,7 +603,7 @@ def dataset_inputs_bootstrap_chi2_central(
     a different value can be specified in the runcard.
     """
     dt, th = dataset_inputs_results
-    diff = np.array(th._rawdata - dt.central_value[:, np.newaxis])
+    diff = np.array(th.rawdata - dt.central_value[:, np.newaxis])
     cchi2 = lambda x, y: calc_chi2(y, x.mean(axis=1))
     chi2_central_resample = bootstrap_values(
         diff,
@@ -727,7 +698,7 @@ def positivity_predictions_data_result(pdf, posdataset):
 def count_negative_points(possets_predictions):
     """Return the number of replicas with negative predictions for each bin
     in the positivity observable."""
-    return np.sum([(r.rawdata < 0).sum(axis=1) for r in possets_predictions], axis=0)
+    return np.sum([(r.rawdata < 0).sum(axis=0) for r in possets_predictions], axis=0)
 
 
 chi2_stat_labels = {
diff --git a/validphys2/src/validphys/tests/baseline/test_dataspecschi2.png b/validphys2/src/validphys/tests/baseline/test_dataspecschi2.png
index 7472b78f63..0115b2daa3 100644
Binary files a/validphys2/src/validphys/tests/baseline/test_dataspecschi2.png and b/validphys2/src/validphys/tests/baseline/test_dataspecschi2.png differ
diff --git a/validphys2/src/validphys/tests/regressions/test_datasetchi2.csv b/validphys2/src/validphys/tests/regressions/test_datasetchi2.csv
index 8d33741347..7f00493b6c 100644
--- a/validphys2/src/validphys/tests/regressions/test_datasetchi2.csv
+++ b/validphys2/src/validphys/tests/regressions/test_datasetchi2.csv
@@ -1,6 +1,6 @@
 		test	test
 		ndata	$\chi^2/ndata$
 group	dataset		
-NMC	NMC	204	1.6064257972017228
-ATLAS	ATLASTTBARTOT	3	1.9383209541765103
-CMS	CMSZDIFF12	28	1.8520436886594904
+NMC	NMC	204	1.6064258130123148
+ATLAS	ATLASTTBARTOT	3	1.9377358549444097
+CMS	CMSZDIFF12	28	1.8519163102572869
diff --git a/validphys2/src/validphys/tests/regressions/test_replicachi2data.csv b/validphys2/src/validphys/tests/regressions/test_replicachi2data.csv
index fed92dac6b..a669de4f65 100644
--- a/validphys2/src/validphys/tests/regressions/test_replicachi2data.csv
+++ b/validphys2/src/validphys/tests/regressions/test_replicachi2data.csv
@@ -1,103 +1,103 @@
 name	Total	NMC	ATLAS	CMS
 npoints	236	204	4	28
-0	1.633164431281718	1.6064257972017228	1.4646799677170539	1.8520436886594904
-1	1.7334088710357112	1.734961856769467	1.1099497616131668	1.8111598477501392
-2	1.68778827233719	1.6453399067644063	2.065152613344062	1.9431457442236326
-3	1.6987178708486776	1.6814301182146243	1.2167258017617955	1.8935275070520494
-4	1.6537880848753	1.620679175623579	1.514372749699407	1.9149266144486803
-5	1.6577037735107656	1.6454811689887392	1.117015026048428	1.8239954275230033
-6	1.6535446769397624	1.6439890548214566	1.1540292226163327	1.7945235601336222
-7	1.690848506785124	1.6949777388063163	1.3206896546250455	1.7136439380821642
-8	1.722111679782665	1.7073333872196168	1.2668130655618413	1.8948247562021334
-9	1.6366119211994845	1.6096259302050373	1.4737139237660677	1.8564952837923716
-10	1.6599397717213789	1.6357991334540212	1.5882755477975263	1.8460593110869656
-11	1.668345767062505	1.6565059983987405	1.7919199070186123	1.73695349019049
-12	1.7178154189496435	1.7074563636898035	1.418145414173511	1.8360985365250682
-13	1.716495702132734	1.71132193676457	1.5322767626714575	1.7805072697381108
-14	1.7425467147007376	1.7329714557322538	1.4174319890603067	1.8587542765626104
-15	1.64044310786435	1.6265082854552777	1.5855998386740953	1.749802995300483
-16	1.661557726572681	1.6587447787248	0.8909779963612721	1.7921348794945864
-17	1.7052372193539282	1.700842932677164	1.723573670346641	1.734633243571396
-18	1.6719119078080142	1.6638605897814052	1.3607457653391302	1.7750238166402903
-19	1.6130225383252468	1.568761881907043	1.8264315671679827	1.9050060309660544
-20	1.8506900570276825	1.8458849948275038	1.2961577816415402	1.9649172638270045
-21	1.6139341914240108	1.583815548848013	1.4818626908545809	1.8522373731304838
-22	1.6996816474924457	1.7204394120543445	0.9511589587802106	1.6553783183575017
-23	1.6928033023915565	1.6761712326724887	1.6229438315790252	1.8239597347465568
-24	1.7131249586938784	1.686905069670985	1.710663225624133	1.904507254870635
-25	1.6368101897760852	1.6263906528587062	1.1344843774701272	1.7844847890749824
-26	1.7761681892360117	1.770018776879609	1.4237754278274368	1.8713128737481715
-27	1.6244279250606297	1.5943841663389933	1.5503462533268624	1.8539012631373752
-28	1.7328041601506325	1.7007419829481452	1.7869796381066878	1.9586606686321744
-29	1.580060371970743	1.5461835710042657	1.2069596503285365	1.8801771678182475
-30	1.6324682267758936	1.612547688513531	1.1424195171903564	1.8476105354853265
-31	1.6076745614174293	1.5880460423569007	1.4760889778463109	1.7694802836542982
-32	1.6206025290769757	1.6050499850907356	1.2982011462916851	1.7799712613746224
-33	1.7735640730921738	1.779356775560228	1.1365669826912466	1.822359682310768
-34	1.597332334515502	1.557894017444925	2.2212178832483884	1.7955421376392933
-35	1.6632195620305688	1.6377982295016988	1.1735879977170232	1.918379493928556
-36	1.6022167994154974	1.5671384218201185	1.4895859069351642	1.873877963679019
-37	1.6223085881457122	1.5990897337236325	1.280468204289206	1.8403088680575101
-38	1.7120969799361094	1.6979197486781215	1.326147968851687	1.8705238092563692
-39	1.6193617601436447	1.5743923499530943	1.8287448078054567	1.9170841704373944
-40	1.7049923143581074	1.6613964971953845	2.381880088571624	1.9259207287988709
-41	1.5938322601370807	1.5724534526025482	0.9703038624537993	1.8386676289862862
-42	1.6196431955840294	1.5919383597818284	1.2681047917285646	1.8717124855508462
-43	1.6735152920612626	1.6423587908608461	1.557506684824832	1.9170853161266457
-44	1.6772042235442266	1.6564533528283856	1.5698541778713182	1.8437248595700546
-45	1.5479597883574368	1.5186200973376842	1.2069809013661184	1.8104316639286775
-46	1.6392941356361497	1.5886581009650398	2.208281316593316	1.9269299338174986
-47	1.7253953793605592	1.7202466235589386	1.5280573974832383	1.7910988833262678
-48	1.7297206748744025	1.6980156851246222	1.7744434332766248	1.9543252061367684
-49	1.568093721907339	1.5232530425694568	1.7658929827122676	1.866533062682633
-50	1.5525538602143887	1.5200417997181	1.6284163368365765	1.7785899471698927
-51	1.6765151911281704	1.6580087457891148	1.7043549850520683	1.8073707508950203
-52	1.6831135357243845	1.669758728082639	1.0365912480541886	1.8727731753528438
-53	1.5010009205633057	1.4336699766955983	1.3116390752971565	2.0186066323517693
-54	1.5425452460692444	1.4879833340395212	1.5852867658031997	1.9339618166095183
-55	1.6962730852353014	1.6654284396198906	1.9025690058549662	1.8915275146304844
-56	1.690682908697377	1.6737669804996216	1.4138561956694147	1.8534742017135901
-57	1.6387017018967793	1.6074557622889447	1.6650462368853587	1.8625871854697766
-58	1.6537047293365406	1.6227941686211176	1.7550825295238084	1.8644277002364416
-59	1.6743055796412083	1.6535803078778097	1.6095502508859874	1.8345547494538572
-60	1.7607179036878984	1.7578708846521396	1.251627048053597	1.854187736039042
-61	1.7112065147096718	1.7094673755827825	1.3069176840456642	1.7816329327290092
-62	1.6760392057425866	1.6445275547369214	1.1999248811986236	1.9736404237187146
-63	1.7648401964721359	1.764419554516598	1.0609413155955352	1.8684618565591433
-64	1.6430979051845722	1.6243586252454583	1.302163581192017	1.8283318481684812
-65	1.526704742353999	1.4950577498604019	1.4384455815477584	1.7698841392082425
-66	1.6223085881457122	1.5990897337236325	1.280468204289206	1.8403088680575101
-67	1.6382019108736001	1.6144363734443945	1.6301782266020446	1.8124970670394653
-68	1.6630312322824974	1.627326621672826	2.057819008210531	1.8667665701632428
-69	1.7235888299300515	1.7108335873411127	1.627832609762473	1.8301993431019734
-70	1.634600771603688	1.6053015596724272	1.4440364915668087	1.8752899271081431
-71	1.6039946631092585	1.5617962838621486	1.8526871064962425	1.8759125057114885
-72	1.608599625679189	1.5908862027577169	1.3561137249874793	1.7737239784915861
-73	1.5918429841463377	1.5497699252028412	1.7656396286339975	1.8735471786650062
-74	1.6668974226763156	1.6476007887206563	1.3512416929491098	1.8525808600285745
-75	1.7739767803927493	1.7741234582593264	1.1282653671703948	1.8651526149680229
-76	1.5341485776112116	1.4801487958594524	2.2995933014729792	1.8182263126794906
-77	1.6207879028642886	1.566560490353848	1.6494267152765514	2.011782077952892
-78	1.6014110063490559	1.5687010433284763	1.2389883627788534	1.8915011145804492
-79	1.618390169594721	1.5985732897790172	1.3178708923949136	1.8057016192805346
-80	1.7039286952694999	1.690631118932874	1.5300461547181594	1.8256514000865363
-81	1.6153472272463718	1.583147892062274	1.0920986319691677	1.9246921829129704
-82	1.5501473720992347	1.5048484827713322	1.420337927786493	1.898726343532917
-83	1.796026945812748	1.8012390234686315	1.3744759444909087	1.8182748087944318
-84	1.6441353118814421	1.6262304929862843	1.5761095481872307	1.7843026729310498
-85	1.7489380132960697	1.7378844464821845	1.569392430447557	1.8551205119184502
-86	1.714996709303176	1.6841598356444387	1.9369792722406798	1.9079535655400477
-87	1.7558507759150124	1.7371573611160664	1.1430508743020027	1.979588498252049
-88	1.6039287515264087	1.56983065566088	1.9306860382050528	1.8056781233068822
-89	1.6419938693980591	1.6367559559810334	1.3686710460232878	1.719201927632786
-90	1.7021197058143331	1.6985942170700976	1.5036495676535802	1.7561582864024405
-91	1.6137527151926292	1.5500116616442616	2.1401226826838564	2.0029561099748467
-92	1.7134425766070225	1.678772855568108	1.4483985152986738	2.003899695791734
-93	1.8242153819893585	1.8257548306283662	1.2598214132543786	1.8936271088673016
-94	1.6333155151589833	1.5883899915574327	1.8319270693757708	1.9322569650821693
-95	1.6178144950080646	1.567217173694048	1.6964396493767955	1.975219956814653
-96	1.5928024780683734	1.5635473746690753	1.1297806308156297	1.8720927810136507
-97	1.6000041690058031	1.5505054784799166	1.5899016684035536	1.9620807000661546
-98	1.6449601173135953	1.6140965055263956	1.9283794416314457	1.8293350997177857
-99	1.6122475097631195	1.576742783494375	1.4402107046238053	1.8955014875981586
-100	1.6457628044745598	1.6294564991751102	1.107293353901103	1.8414900931667597
+0	1.6331414605961554	1.6064258130123148	1.4642155397440857	1.851916310257288
+1	1.7334090365532264	1.734962043323776	1.109948798376749	1.8111600212501475
+2	1.6877884536672854	1.6453401923939792	2.065149379329199	1.9431456535639562
+3	1.6987177907128304	1.6814300483245213	1.2167235538564396	1.8935276619499943
+4	1.6537882304597358	1.6206794198129082	1.5143707501014874	1.9149263480806569
+5	1.657704094906477	1.6454815520518558	1.1170136902028251	1.8239955363763802
+6	1.6535444344982564	1.6439888069898139	1.1540274992869584	1.7945235685185232
+7	1.6908488964010466	1.6949782613861923	1.3206871997065657	1.7136437653227667
+8	1.7221117793168024	1.7073335549373894	1.2668104711359693	1.8948247438212158
+9	1.6366119834769124	1.6096260453398075	1.473712423033797	1.856495184253406
+10	1.659940184824683	1.635799670902719	1.588272788114648	1.8460592715004263
+11	1.6683459190719685	1.6565061881106242	1.7919180189232864	1.7369536589544328
+12	1.7178157546717383	1.7074567962836644	1.4181426918961695	1.836098603324215
+13	1.716495766243632	1.7113220642754658	1.5322745582378068	1.7805071960125347
+14	1.7425466646065497	1.7329714408352708	1.4174294913223486	1.8587543196950382
+15	1.6404432869597365	1.6265085381351183	1.5855982433010838	1.7498028917760478
+16	1.6615573279924452	1.6587443746631547	0.8909766503971864	1.7921346561908817
+17	1.7052376891371723	1.7008434882599703	1.723571900842281	1.7346334081417716
+18	1.6719120691278573	1.6638607953367437	1.3607436977818923	1.7750239740839657
+19	1.613022802794292	1.5687621999135686	1.8264293279129227	1.905006263051188
+20	1.8506898474160518	1.8458847952507411	1.2961558920638951	1.9649172210993366
+21	1.6139342860350945	1.5838156819370806	1.4818614482670958	1.8522373784303365
+22	1.6996814373159437	1.7204391927861087	0.9511578379601898	1.655378304512707
+23	1.692803617766823	1.676171636000037	1.6229418116006402	1.8239597429485772
+24	1.7131251359911894	1.686905321858746	1.7106610295319293	1.9045072255931714
+25	1.6368102934945292	1.6263908361520953	1.1344820679308685	1.784484657784212
+26	1.7761678251548632	1.7700184404698238	1.4237723277616203	1.8713126989163273
+27	1.6244278267962076	1.5943841151893652	1.5503438175289963	1.8539011555413756
+28	1.7328040172990973	1.700741869031167	1.7869767063887179	1.9586607133812157
+29	1.5800604839603773	1.5461837217018417	1.2069577326128216	1.880177287750787
+30	1.6324681919097153	1.612547707995762	1.1424179131898138	1.8476103288142185
+31	1.607674369706397	1.5880458815229432	1.4760862036609292	1.7694802359066302
+32	1.6206027204928388	1.605050271148119	1.2981988532110809	1.7799711181874804
+33	1.773563907994167	1.779356634055176	1.1365643674204111	1.8223596953459236
+34	1.5973323690956018	1.5578941707519829	2.2212137553099227	1.7955419018542085
+35	1.6632202782908945	1.637799088179375	1.1735866589147328	1.9183794661571316
+36	1.6022165691797223	1.5671381536485012	1.4895846454956718	1.8738781571477698
+37	1.6223084864434523	1.5990896439417261	1.2804662364817754	1.8403089460934123
+38	1.7120967709150565	1.6979195269489693	1.326145536357827	1.870524010461868
+39	1.6193616450506452	1.57439224522193	1.8287421073169994	1.9170843491932335
+40	1.7049922104767428	1.6613964627270275	2.381877322608255	1.925920499491593
+41	1.5938325521064143	1.572453807552316	0.970302242252107	1.8386677352654597
+42	1.6196429364071514	1.5919380830753238	1.2681034146992118	1.8717125137830297
+43	1.6735152189707765	1.6423587572506044	1.557504332292557	1.917085281028917
+44	1.677203841320786	1.6564529677486477	1.5698518477272483	1.843724776431155
+45	1.5479594468185343	1.5186197513713733	1.2069787364086781	1.8104316151349715
+46	1.6392941239480796	1.5886581532141189	2.208278196687121	1.926929900332788
+47	1.7253951918147818	1.7202464401597746	1.5280554477890402	1.7910989173049416
+48	1.7297207556439842	1.6980158287288005	1.7744412647645693	1.9543251504373793
+49	1.5680938248879044	1.523253230985689	1.7658909257774207	1.8665328517626851
+50	1.5525536064711594	1.5200415339606121	1.6284144798655795	1.778590009991659
+51	1.676515550905596	1.6580091778992072	1.7043532958734795	1.8073708763853022
+52	1.6831136663834798	1.6697588923581719	1.0365893774337271	1.872773346989258
+53	1.5010009180846036	1.4336700079261697	1.3116366641605577	2.0186067283709135
+54	1.542545356916207	1.4879834892236081	1.5852844792578435	1.9339619469134803
+55	1.6962726312543277	1.665427945689167	1.902566187404671	1.8915276894933064
+56	1.6906827163282356	1.6737668230523257	1.4138531349916628	1.8534741646722332
+57	1.6387017064660885	1.6074558030941164	1.6650435522964906	1.862587310200399
+58	1.6537049913056714	1.6227945488852065	1.7550794096713644	1.8644275834596755
+59	1.6743054819772532	1.653580252097564	1.6095477186920686	1.8345546944271567
+60	1.760717962262651	1.7578709735097442	1.2516251105595473	1.8541878591342729
+61	1.7112061457290064	1.70946698914875	1.3069157010114216	1.7816329214876736
+62	1.6760395309755631	1.6445279689735608	1.199922754227893	1.9736404508112462
+63	1.7648405518025423	1.764419975860034	1.0609399045624792	1.8684619832751128
+64	1.643097893117864	1.6243586788019808	1.3021610081017991	1.8283317238501635
+65	1.5267050565132614	1.495058150632585	1.438443534733451	1.7698841596124473
+66	1.6223084864434523	1.5990896439417261	1.2804662364817754	1.8403089460934123
+67	1.6382020526715084	1.6144365545196706	1.6301756722086964	1.8124973078438722
+68	1.6630309961175795	1.6273264451809057	2.0578155613271347	1.8667663579119824
+69	1.7235888968944642	1.7108337002459408	1.6278304630037084	1.8301993916038146
+70	1.6346012094259845	1.6053021309289885	1.4440341757319628	1.8752897861461022
+71	1.6039943852756267	1.5617960046403603	1.8526844960263371	1.8759125712253213
+72	1.6085997177416667	1.5908863719277442	1.356111823911515	1.7737237935045518
+73	1.5918426078675785	1.549769525608617	1.765636357099144	1.873547385864074
+74	1.6668975328016988	1.6476009491016603	1.3512393598159287	1.8525809530428037
+75	1.7739771273858427	1.774123891433172	1.1282642608772993	1.8651525416850916
+76	1.534148507374288	1.480148785522714	2.299589358068434	1.818226359336592
+77	1.6207881011370513	1.5665607783792113	1.6494241205241973	2.0117820213174356
+78	1.6014112826764793	1.5687014104441326	1.2389866663656035	1.8915010112708424
+79	1.6183903753681967	1.5985735939420136	1.31786788462861	1.8057015672931855
+80	1.7039281449930659	1.6906305465912053	1.530042988268256	1.825651384310167
+81	1.6153470453271934	1.5831477255608042	1.0920969315185147	1.9246921055978423
+82	1.550147327874484	1.5048484618542393	1.4203353054935601	1.898726497790684
+83	1.796026330456324	1.8012384085493949	1.374472472263034	1.8182745983772783
+84	1.6441355146605652	1.6262307965013056	1.5761068082679037	1.7843025621626933
+85	1.7489385095427001	1.7378851115215144	1.569388920861313	1.8551203506515386
+86	1.7149965294782474	1.684159684138558	1.9369763955439816	1.9079535646580241
+87	1.7558511808377928	1.7371578723123293	1.1430486440552472	1.9795885053493925
+88	1.60392879188972	1.5698307808261078	1.9306828018708304	1.8056780139273083
+89	1.6419937788873573	1.6367559127549003	1.3686691251536447	1.7192017541000761
+90	1.7021195976333319	1.6985941849045854	1.5036475642040819	1.7561578951469468
+91	1.6137530834181903	1.5500121223068863	2.1401203891314338	2.002956184984369
+92	1.7134421727269178	1.6787724457282462	1.4483964149900914	2.003899577679644
+93	1.8242157601462674	1.8257553166701763	1.2598197029719003	1.8936269993541253
+94	1.6333154905545224	1.5883899834841848	1.831924912678808	1.9322571246206564
+95	1.6178142924268402	1.5672169871582218	1.6964380697827446	1.9752198340473586
+96	1.5928022801018331	1.5635472024642867	1.129778609475359	1.8720926558363122
+97	1.6000044124178256	1.5505058092904849	1.5898997018858283	1.9620806224215934
+98	1.6449604978934307	1.6140969815229855	1.9283777115501315	1.8293350866414326
+99	1.6122478293595508	1.576743205701731	1.4402089041553603	1.8955013624671206
+100	1.6457631138326914	1.629456917858709	1.1072910758816195	1.8414899756361456
diff --git a/validphys2/src/validphys/tests/test_closuretest.py b/validphys2/src/validphys/tests/test_closuretest.py
index a55eb4063a..68fc91dc6e 100644
--- a/validphys2/src/validphys/tests/test_closuretest.py
+++ b/validphys2/src/validphys/tests/test_closuretest.py
@@ -13,7 +13,7 @@ class TestResult:
     """class for testing base level estimators which expect a results object"""
     def __init__(self, central_value, rawdata=None):
         self.central_value = central_value
-        self._rawdata = rawdata
+        self.rawdata = rawdata
         self.ndata = len(central_value)
         self.sqrtcovmat = np.identity(self.ndata)
 
diff --git a/validphys2/src/validphys/tests/test_pyfkdata.py b/validphys2/src/validphys/tests/test_pyfkdata.py
index 0f6b90b6a2..37c45720c9 100644
--- a/validphys2/src/validphys/tests/test_pyfkdata.py
+++ b/validphys2/src/validphys/tests/test_pyfkdata.py
@@ -68,7 +68,7 @@ def test_predictions(pdf_name):
         ds = l.check_dataset(**daset, theoryid=THEORYID)
         preds = predictions(ds, pdf)
         core_predictions = ThPredictionsResult.from_convolution(pdf, ds)
-        assert_allclose(preds.values, core_predictions._rawdata, atol=1e-8, rtol=1e-3)
+        assert_allclose(preds.values, core_predictions.rawdata, atol=1e-8, rtol=1e-3)
         # Now check that the stats class does the right thing with the data
         cv_predictions = central_predictions(ds, pdf).squeeze()
         stats_predictions = pdf.stats_class(preds.T)
@@ -90,12 +90,12 @@ def test_positivity(pdf_name):
         ps = l.check_posset(setname=posset, theoryID=THEORYID, postlambda=1e6)
         preds = predictions(ps, pdf)
         core_predictions = PositivityResult.from_convolution(pdf, ps)
-        assert_allclose(preds.values, core_predictions.rawdata.T)
+        assert_allclose(preds.values, core_predictions.rawdata)
         # Now do the same with the API
         api_predictions = API.positivity_predictions_data_result(
             theoryid=THEORYID, pdf=pdf_name, posdataset={"dataset": posset, "maxlambda": 1e6}
         )
-        assert_allclose(preds.values, api_predictions.rawdata.T)
+        assert_allclose(preds.values, api_predictions.rawdata)
         # And now check that the results are correct for any kind of PDF
         cv_predictions = central_predictions(ps, pdf).squeeze()
         assert_allclose(cv_predictions, api_predictions.central_value, atol=1e-3)